Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F86077157
FE_Engine.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Fri, Oct 4, 03:07
Size
97 KB
Mime Type
text/x-c
Expires
Sun, Oct 6, 03:07 (2 d)
Engine
blob
Format
Raw Data
Handle
21341961
Attached To
rLAMMPS lammps
FE_Engine.cpp
View Options
#include "FE_Engine.h"
#include "ATC_Transfer.h"
#include "FE_Element.h"
#include "Function.h"
#include "PhysicsModel.h"
#include "KernelFunction.h"
#include "Utility.h"
#include <stdio.h>
#include <stdlib.h>
#include <map>
#include <string>
using namespace std;
namespace ATC{
static const double tol_sparse = 1.e-30;//tolerance for compaction from dense
//-----------------------------------------------------------------
FE_Engine::FE_Engine(MPI_Comm comm)
: communicator_(comm),
feMesh_(NULL),
initialized_(false),
outputManager_()
{
// Nothing to do here
}
//-----------------------------------------------------------------
FE_Engine::~FE_Engine()
{
if (feMesh_) delete feMesh_;
}
//-----------------------------------------------------------------
void FE_Engine::initialize()
{
if (!feMesh_) throw ATC_Error("FE_Engine has no mesh");
if (!initialized_) {
// set up work spaces
nNodesPerElement_ = feMesh_->num_nodes_per_element();
nIPsPerElement_ = feMesh_->num_ips_per_element();
nIPsPerFace_ = feMesh_->num_ips_per_face();
nSD_ = feMesh_->num_spatial_dimensions();
nElems_ = feMesh_->num_elements();
nNodesUnique_ = feMesh_->num_nodes_unique();
nNodes_ = feMesh_->num_nodes();
// arrays & matrices
_weights_.reset(nIPsPerElement_,nIPsPerElement_);
_N_.reset(nIPsPerElement_,nNodesPerElement_);
_dN_.assign(nSD_, DENS_MAT(nIPsPerElement_,nNodesPerElement_));
_Nw_.reset(nIPsPerElement_,nNodesPerElement_);
_dNw_.assign(nSD_, DENS_MAT(nIPsPerElement_,nNodesPerElement_));
_Bfluxes_.assign(nSD_, DENS_MAT());
// faces
_fweights_.reset(nIPsPerElement_,nIPsPerElement_);
_fN_.reset(nIPsPerFace_,nNodesPerElement_);
_fdN_.assign(nSD_, DENS_MAT(nIPsPerFace_, nNodesPerElement_));
_nN_.assign(nSD_, DENS_MAT(nIPsPerFace_, nNodesPerElement_));
// remove specified elements
if (nullElements_.size() > 0) delete_elements(nullElements_);
initialized_ = true;
}
}
void FE_Engine::partition_mesh()
{
if (is_partitioned()) return;
feMesh_->partition_mesh();
// now do all FE_Engine data structure partitioning
// partition nullElements_
/*for (vector<int>::iterator elemsIter = feMesh_->processor_elts().begin();
elemsIter != feMesh_->processor_elts().end();
++elemsIter)
{
if (nullElements_.find(*elemsIter) != nullElements_.end()) {
myNullElements_.insert(map_elem_to_myElem(*elemsIter));
}
}*/
}
void FE_Engine::departition_mesh()
{
if (!is_partitioned()) return;
feMesh_->departition_mesh();
}
void FE_Engine::set_quadrature(FeIntQuadrature type, bool temp) const
{
if (!feMesh_) throw ATC_Error("FE_Engine has no mesh");
feMesh_->set_quadrature(type);
if (!temp) quadrature_ = type;
int nIPsPerElement_new = feMesh_->num_ips_per_element();
int nIPsPerFace_new = feMesh_->num_ips_per_face();
if (nIPsPerElement_ != nIPsPerElement_new) {
// arrays & matrices
nIPsPerElement_ = nIPsPerElement_new;
_weights_.resize(nIPsPerElement_,nIPsPerElement_);
_N_.resize(nIPsPerElement_,nNodesPerElement_);
_dN_.assign(nSD_, DENS_MAT(nIPsPerElement_,nNodesPerElement_));
_Nw_.reset(nIPsPerElement_,nNodesPerElement_);
_dNw_.assign(nSD_, DENS_MAT(nIPsPerElement_,nNodesPerElement_));
}
if (nIPsPerFace_ != nIPsPerFace_new) {
// faces
nIPsPerFace_ = nIPsPerFace_new;
_fweights_.reset(nIPsPerElement_,nIPsPerElement_);
_fN_.reset(nIPsPerFace_,nNodesPerElement_);
_fdN_.assign(nSD_, DENS_MAT(nIPsPerFace_, nNodesPerElement_));
_nN_.assign(nSD_, DENS_MAT(nIPsPerFace_, nNodesPerElement_));
}
}
//-----------------------------------------------------------------
bool FE_Engine::modify(int narg, char **arg)
{
bool match = false;
/*! \page man_mesh fix_modify AtC mesh create
\section syntax
fix_modify AtC mesh create <nx> <ny> <nz> <region-id>
<f|p> <f|p> <f|p>
- nx ny nz = number of elements in x, y, z
- region-id = id of region that is to be meshed
- f p p = perioidicity flags for x, y, z
fix_modify AtC mesh quadrature <quad> where
- quad = one of <nodal|gauss1|gauss2|gauss3|face>
- when a mesh is created it defaults to gauss2, use this call
to change it after the fact
\section examples
<TT> fix_modify AtC create mesh 10 1 1 feRegion p p p </TT>
\section description
Creates a uniform mesh in a rectangular region
\section restrictions
creates only uniform rectangular grids in a rectangular region
\section related
\section default
none
*/
int argIdx = 0;
if (strcmp(arg[argIdx],"mesh")==0) {
argIdx++;
// create mesh
if (strcmp(arg[argIdx],"create")==0) {
if (feMesh_) throw ATC_Error("FE Engine already has a mesh");
argIdx++;
int nx = atoi(arg[argIdx++]);
int ny = atoi(arg[argIdx++]);
int nz = atoi(arg[argIdx++]);
string box = arg[argIdx++];
Array<bool> periodicity(3);
periodicity(0) = (strcmp(arg[argIdx++],"p")==0) ? true : false;
periodicity(1) = (strcmp(arg[argIdx++],"p")==0) ? true : false;
periodicity(2) = (strcmp(arg[argIdx++],"p")==0) ? true : false;
if (argIdx < narg ) {
Array<double> dx(nx),dy(ny),dz(nz);
dx = 0;
dy = 0;
dz = 0;
double x[3] = {0,0,0};
while (argIdx < narg) {
if (strcmp(arg[argIdx++],"dx")==0) {
// parse relative values for each element
if (is_numeric(arg[argIdx])) {
for (int i = 0; i < nx; ++i) {
if (is_numeric(arg[argIdx])) { dx(i) = atof(arg[argIdx++]); }
else throw ATC_Error("not enough element partitions");
}
}
// construct relative values from a density function
// evaluate for a domain (0,1)
else {
XT_Function * f = XT_Function_Mgr::instance()->function(&(arg[argIdx]),narg-argIdx);
for (int i = 0; i < nx; ++i) {
x[0] = (i+0.5)/nx; dx(i) = f->f(x,0.);
}
}
break ;
}
else if (strcmp(arg[argIdx++],"dy")==0) {
for (int i = 0; i < ny; ++i) { dy(i) = atof(arg[argIdx++]); }
}
else if (strcmp(arg[argIdx++],"dz")==0) {
for (int i = 0; i < nz; ++i) { dz(i) = atof(arg[argIdx++]); }
}
}
create_mesh(dx, dy, dz, box.c_str(), periodicity);
}
else {
create_mesh(nx, ny, nz, box.c_str(), periodicity);
}
quadrature_ = GAUSS2;
match = true;
}
else if (strcmp(arg[argIdx],"quadrature")==0) {
argIdx++;
string quadStr = arg[argIdx];
FeIntQuadrature quadEnum = string_to_FIQ(quadStr);
set_quadrature(quadEnum);
match = true;
}
/** Example command for reading a mesh from a file:
fix_modify atc mesh read fileName */
// read mesh
else if (strcmp(arg[argIdx],"read")==0) {
argIdx++;
string meshFile = arg[argIdx++];
Array<bool> periodicity(3);
periodicity = false;
if (argIdx < narg) {
periodicity(0) = (strcmp(arg[argIdx++],"p")==0) ? true : false;
periodicity(1) = (strcmp(arg[argIdx++],"p")==0) ? true : false;
periodicity(2) = (strcmp(arg[argIdx++],"p")==0) ? true : false;
}
read_mesh(meshFile,periodicity);
if (periodicity(0) || periodicity(1) || periodicity(2)) {
meshFile = "periodic_"+meshFile;
stringstream ss;
ss << "writing periodicity corrected mesh: " << meshFile;
print_msg(communicator_,ss.str());
feMesh_->write_mesh(meshFile);
feMesh_->output(meshFile);
}
match = true;
}
/** Example command for reading a mesh from a file:
fix_modify atc mesh write fileName */
// write mesh
else if (strcmp(arg[argIdx],"write")==0) {
argIdx++;
string meshFile = arg[argIdx];
feMesh_->write_mesh(meshFile);
match = true;
}
/*! \page man_mesh_delete_elements fix_modify AtC mesh delete_elements
\section syntax
fix_modify AtC mesh delete_elements <element_set>
- <element_set> = name of an element set
\section examples
<TT> fix_modify AtC delete_elements gap </TT>
\section description
Deletes a group of elements from the mesh.
\section restrictions
\section related
\section default
none
*/
else if (strcmp(arg[argIdx],"delete_elements")==0) {
argIdx++;
string esetName = arg[argIdx];
set<int> elemSet = feMesh_->elementset(esetName);
nullElements_.insert(elemSet.begin(), elemSet.end());
match = true;
}
else if (strcmp(arg[argIdx],"cut")==0) {
argIdx++;
string fsetName = arg[argIdx++];
set<PAIR> faceSet = feMesh_->faceset(fsetName);
cutFaces_.insert(faceSet.begin(), faceSet.end());
if (narg > argIdx && strcmp(arg[argIdx],"edge")==0) {
argIdx++;
string nsetName = arg[argIdx];
set<int> nodeSet = feMesh_->nodeset(nsetName);
cutEdge_.insert(nodeSet.begin(), nodeSet.end());
}
// cut mesh
if (cutFaces_.size() > 0) cut_mesh(cutFaces_,cutEdge_);
match = true;
}
else if (strcmp(arg[argIdx],"lammps_partition")==0) {
feMesh_->set_lammps_partition(true);
match = true;
}
else if (strcmp(arg[argIdx],"data_decomposition")==0) {
feMesh_->set_data_decomposition(true);
match = true;
}
else {
if ( ! feMesh_ ) throw ATC_Error("need mesh for parsing");
match = feMesh_->modify(narg,arg);
}
}
// FE_Mesh
else {
if ( ! feMesh_ ) throw ATC_Error("need mesh for parsing");
match = feMesh_->modify(narg,arg);
}
return match;
}
//-----------------------------------------------------------------
void FE_Engine::finish()
{
// Nothing to do
}
//-----------------------------------------------------------------
// write geometry
//-----------------------------------------------------------------
void FE_Engine::initialize_output(int rank,
string outputPrefix, set<int> otypes)
{
outputManager_.initialize(outputPrefix, otypes);
if (!feMesh_) throw ATC_Error("output needs mesh");
if (!initialized_) initialize();
if (!feMesh_->coordinates() || !feMesh_->connectivity())
throw ATC_Error("output mesh not properly initialized");
if (!feMesh_->coordinates()->nCols() ||
!feMesh_->connectivity()->nCols())
throw ATC_Error("output mesh is empty");
if (rank == 0)
outputManager_.write_geometry(feMesh_->coordinates(),
feMesh_->connectivity());
outputManager_.print_custom_names();
}
//-----------------------------------------------------------------
// write geometry
//-----------------------------------------------------------------
void FE_Engine::write_geometry(void)
{
outputManager_.write_geometry(feMesh_->coordinates(),
feMesh_->connectivity());
}
// -------------------------------------------------------------
// write data
// -------------------------------------------------------------
void FE_Engine::write_data(double time,
FIELDS &soln,
OUTPUT_LIST *data)
{
outputManager_.write_data(
time, &soln, data,
(feMesh_->node_map())->data());
}
// -------------------------------------------------------------
// write data
// -------------------------------------------------------------
void FE_Engine::write_data(double time, OUTPUT_LIST *data)
{
outputManager_.write_data(
time, data,
feMesh_->node_map()->data());
}
// -------------------------------------------------------------
// amend mesh for deleted elements
// -------------------------------------------------------------
void FE_Engine::delete_elements(const set<int> & elementList)
{
feMesh_->delete_elements(elementList);
}
// -------------------------------------------------------------
// amend mesh for cut at specified faces
// -------------------------------------------------------------
void FE_Engine::cut_mesh(const set<PAIR> &faceSet,
const set<int> &nodeSet)
{
feMesh_->cut_mesh(faceSet,nodeSet);
}
// -------------------------------------------------------------
// interpolate one value
// -------------------------------------------------------------
DENS_VEC FE_Engine::interpolate_field(const DENS_VEC & x, const FIELD & f) const
{
DENS_VEC N;
Array<int> nodelist;
feMesh_->shape_functions(x, N, nodelist);
const DENS_MAT &vI = f.quantity();
int dof = vI.nCols();
DENS_MAT vIe(nNodesPerElement_, dof);
for (int i = 0; i < nNodesPerElement_; i++)
for (int j = 0; j < dof; j++)
vIe(i,j) = vI(nodelist(i),j);
DENS_VEC vP;
vP = N*vIe;
return vP;
}
// -------------------------------------------------------------
// interpolate fields and gradients
// Currently, this function will break if called with an unowned ielem.
// Currently, this function is only called with owned ielems.
// -------------------------------------------------------------
void FE_Engine::interpolate_fields(
const int ielem,
const FIELDS &fields,
AliasArray<int> & conn,
DENS_MAT & N,
DENS_MAT_VEC & dN,
DIAG_MAT & _weights_,
FIELD_MATS &fieldsAtIPs,
GRAD_FIELD_MATS &gradFieldsAtIPs) const
{
// evaluate shape functions
feMesh_->shape_function(ielem, N, dN, _weights_);
// get connectivity
_conn_ = feMesh_->element_connectivity_unique(ielem);
// compute fields and gradients of fields ips of this element
FIELD_MATS localElementFields;
for (_fieldItr_ = fields.begin(); _fieldItr_ != fields.end(); _fieldItr_++) {
// field values at all nodes
_fieldName_ = _fieldItr_->first;
const DENS_MAT &vI = (_fieldItr_->second).quantity();
int dof = vI.nCols();
// field values at integration points -> to be computed
DENS_MAT &vP = fieldsAtIPs[_fieldName_];
// gradients of field at integration points -> to be computed
DENS_MAT_VEC &dvP = gradFieldsAtIPs[_fieldName_];
if (_fieldName_ == ELECTRON_WAVEFUNCTION_ENERGIES ) {
vP = vI;
continue;
}
// field values at element nodes
DENS_MAT &vIe = localElementFields[_fieldName_];
// gather local field
vIe.reset(nNodesPerElement_, dof);
for (int i = 0; i < nNodesPerElement_; i++)
for (int j = 0; j < dof; j++)
vIe(i,j) = vI(conn(i),j);
// interpolate field at integration points
vP = N*vIe;
// gradients
dvP.assign(nSD_, DENS_MAT(nIPsPerElement_, dof));
for (int j = 0; j < nSD_; ++j) dvP[j] = dN[j]*vIe;
}
}
// -------------------------------------------------------------
// interpolate fields
// Currently, this function will break if called with an unowned ielem.
// Currently, this function is only called with owned ielems.
// -------------------------------------------------------------
void FE_Engine::interpolate_fields(
const int ielem,
const FIELDS &fields,
AliasArray<int> & conn,
DENS_MAT & N,
DIAG_MAT & _weights_,
FIELD_MATS &fieldsAtIPs) const
{
// evaluate shape functions
feMesh_->shape_function(ielem, N, _weights_);
// get connectivity
_conn_ = feMesh_->element_connectivity_unique(ielem);
// compute fields and gradients of fields ips of this element
FIELD_MATS localElementFields;
for (_fieldItr_ = fields.begin(); _fieldItr_ != fields.end(); _fieldItr_++) {
// field values at all nodes
_fieldName_ = _fieldItr_->first;
const DENS_MAT &vI = (_fieldItr_->second).quantity();
int dof = vI.nCols();
// field values at integration points -> to be computed
DENS_MAT &vP = fieldsAtIPs[_fieldName_];
// field values at element nodes
DENS_MAT &vIe = localElementFields[_fieldName_];
if (_fieldName_ == ELECTRON_WAVEFUNCTION_ENERGIES ) {
vP = vI;
continue;
}
// gather local field
vIe.reset(nNodesPerElement_, dof);
for (int i = 0; i < nNodesPerElement_; i++)
for (int j = 0; j < dof; j++)
vIe(i,j) = vI(conn(i),j);
// interpolate field at integration points
vP = N*vIe;
}
}
// -------------------------------------------------------------
// compute dimensionless stiffness matrix using native quadrature
// -------------------------------------------------------------
void FE_Engine::stiffness_matrix(SPAR_MAT &matrix) const
{
// assemble consistent mass
matrix.reset(nNodesUnique_,nNodesUnique_);// zero since partial fill
DENS_MAT elementMassMatrix(nNodesPerElement_,nNodesPerElement_);
vector<int> myElems = feMesh_->owned_elts();
for (vector<int>::iterator elemsIter = myElems.begin();
elemsIter != myElems.end();
++elemsIter)
{
int ielem = *elemsIter;
// evaluate shape functions
feMesh_->shape_function(ielem, _N_, _dN_, _weights_); // _N_ unused
// perform quadrature
elementMassMatrix = _dN_[0].transMat(_weights_*_dN_[0]);
for (int i = 1; i < nSD_; ++i) {
elementMassMatrix += _dN_[i].transMat(_weights_*_dN_[i]);
}
// get connectivity
_conn_ = feMesh_->element_connectivity_unique(ielem);
for (int i = 0; i < nNodesPerElement_; ++i)
{
int inode = _conn_(i);
for (int j = 0; j < nNodesPerElement_; ++j)
{
int jnode = _conn_(j);
matrix.add(inode, jnode, elementMassMatrix(i,j));
}
}
}
#ifdef ISOLATE_FE
sparse_allsum(communicator_,matrix);
#else
LammpsInterface::instance()->sparse_allsum(matrix);
#endif
matrix.compress();
}
// -------------------------------------------------------------
// compute tangent using native quadrature for one (field,field) pair
// -------------------------------------------------------------
void FE_Engine::compute_tangent_matrix(
const Array2D<bool> &rhsMask,
const pair<FieldName,FieldName> row_col,
const FIELDS &fields,
const PhysicsModel * physicsModel,
const Array<int> & elementMaterials,
SPAR_MAT &tangent,
const DenseMatrix<bool> *elementMask) const
{
tangent.reset(nNodesUnique_,nNodesUnique_);
FieldName rowField = row_col.first;
FieldName colField = row_col.second;
bool BB = rhsMask(rowField,FLUX);
bool NN = rhsMask(rowField,SOURCE);
DENS_MAT elementMatrix(nNodesPerElement_,nNodesPerElement_);
DENS_MAT coefsAtIPs;
vector<int> myElems = feMesh_->owned_elts();
for (vector<int>::iterator elemsIter = myElems.begin();
elemsIter != myElems.end();
++elemsIter)
{
int ielem = *elemsIter;
// if element is masked, skip it
if (elementMask && !(*elementMask)(ielem,0)) continue;
// material id
int imat = elementMaterials(ielem);
const Material * mat = physicsModel->material(imat);
// interpolate fields and gradients (nonlinear only)
interpolate_fields(ielem,fields,_conn_,_N_,_dN_,_weights_,
_fieldsAtIPs_,_gradFieldsAtIPs_);
// evaluate Physics model
if (! (physicsModel->null(rowField,imat)) ) {
if (BB && physicsModel->weak_equation(rowField)->
has_BB_tangent_coefficients() ) {
physicsModel->weak_equation(rowField)->
BB_tangent_coefficients(colField, _fieldsAtIPs_, mat, coefsAtIPs);
DIAG_MAT D(column(coefsAtIPs,0));
D = _weights_*D;
elementMatrix = _dN_[0].transMat(D*_dN_[0]);
for (int i = 1; i < nSD_; i++) {
elementMatrix += _dN_[i].transMat(D*_dN_[i]);
}
}
else {
elementMatrix.reset(nNodesPerElement_,nNodesPerElement_);
}
if (NN && physicsModel->weak_equation(rowField)->
has_NN_tangent_coefficients() ) {
physicsModel->weak_equation(rowField)->
NN_tangent_coefficients(colField, _fieldsAtIPs_, mat, coefsAtIPs);
DIAG_MAT D(column(coefsAtIPs,0));
D = _weights_*D;
elementMatrix += _N_.transMat(D*_N_);
}
// assemble
for (int i = 0; i < nNodesPerElement_; ++i)
{
int inode = _conn_(i);
for (int j = 0; j < nNodesPerElement_; ++j)
{
int jnode = _conn_(j);
tangent.add(inode, jnode, elementMatrix(i,j));
}
}
}
}
#ifdef ISOLATE_FE
sparse_allsum(communicator_,tangent);
#else
LammpsInterface::instance()->sparse_allsum(tangent);
#endif
tangent.compress();
}
// -------------------------------------------------------------
// compute tangent using given quadrature for one (field,field) pair
// -------------------------------------------------------------
void FE_Engine::compute_tangent_matrix(const Array2D<bool> &rhsMask,
const pair<FieldName,FieldName> row_col,
const FIELDS &fields,
const PhysicsModel * physicsModel,
const Array<set<int> > & pointMaterialGroups,
const DIAG_MAT &weights,
const SPAR_MAT &N,
const SPAR_MAT_VEC &dN,
SPAR_MAT &tangent,
const DenseMatrix<bool> *elementMask ) const
{
int nn = nNodesUnique_;
FieldName rowField = row_col.first;
FieldName colField = row_col.second;
bool BB = rhsMask(rowField,FLUX);
bool NN = rhsMask(rowField,SOURCE);
DENS_MAT K(nn,nn);
DENS_MAT coefsAtIPs;
int nips = weights.nCols();
if (nips>0) {
// compute fields and gradients of fields at given ips
GRAD_FIELD_MATS gradFieldsAtIPs;
FIELD_MATS fieldsAtIPs;
for (_fieldItr_ = fields.begin(); _fieldItr_ != fields.end(); _fieldItr_++) {
_fieldName_ = _fieldItr_->first;
const DENS_MAT & field = (_fieldItr_->second).quantity();
int dof = field.nCols();
gradFieldsAtIPs[_fieldName_].assign(nSD_,DENS_MAT(nips,dof));
fieldsAtIPs[_fieldName_] = N*field;
for (int j = 0; j < nSD_; ++j) {
gradFieldsAtIPs[_fieldName_][j] = (*dN[j])*field;
}
}
// treat single material point sets specially
int nMatls = pointMaterialGroups.size();
int atomMatls = 0;
for (int imat = 0; imat < nMatls; imat++) {
const set<int> & indices = pointMaterialGroups(imat);
if ( indices.size() > 0) atomMatls++;
}
bool singleMaterial = ( atomMatls == 1 );
if (! singleMaterial ) throw ATC_Error("FE_Engine::compute_tangent_matrix-given quadrature can not handle multiple atom material currently");
if (singleMaterial)
{
int imat = 0;
const Material * mat = physicsModel->material(imat);
// evaluate Physics model
if (! (physicsModel->null(rowField,imat)) ) {
if (BB && physicsModel->weak_equation(rowField)->
has_BB_tangent_coefficients() ) {
physicsModel->weak_equation(rowField)->
BB_tangent_coefficients(colField, fieldsAtIPs, mat, coefsAtIPs);
DIAG_MAT D(column(coefsAtIPs,0));
D = weights*D;
K = (*dN[0]).transMat(D*(*dN[0]));
for (int i = 1; i < nSD_; i++) {
K += (*dN[i]).transMat(D*(*dN[i]));
}
}
if (NN && physicsModel->weak_equation(rowField)->
has_NN_tangent_coefficients() ) {
physicsModel->weak_equation(rowField)->
NN_tangent_coefficients(colField, fieldsAtIPs, mat, coefsAtIPs);
DIAG_MAT D(column(coefsAtIPs,0));
D = weights*D;
K += N.transMat(D*N);
}
}
}
}
// share information between processors
int count = nn*nn;
DENS_MAT K_sum(nn,nn);
allsum(communicator_, K.ptr(), K_sum.ptr(), count);
// create sparse from dense
tangent.reset(K_sum,tol_sparse);
tangent.compress();
}
// -------------------------------------------------------------
// compute a consistent mass matrix for native quadrature
// -------------------------------------------------------------
void FE_Engine::compute_mass_matrix(
const Array<FieldName>& field_mask,
const FIELDS &fields,
const PhysicsModel * physicsModel,
const Array<int> & elementMaterials,
CON_MASS_MATS & massMatrices,
const DenseMatrix<bool> *elementMask) const
{
int nfields = field_mask.size();
vector<FieldName> usedFields;
DENS_MAT elementMassMatrix(nNodesPerElement_,nNodesPerElement_);
// (JAT, 04/21/11) FIX THIS
DENS_MAT capacity;
// zero, use incoming matrix as template if possible
for (int j = 0; j < nfields; ++j) {
_fieldName_ = field_mask(j);
SPAR_MAT & M = massMatrices[_fieldName_].set_quantity();
// compresses 2May11
if (M.has_template()) { M = 0; }
else { M.reset(nNodesUnique_,nNodesUnique_); }
M.reset(nNodesUnique_,nNodesUnique_);
}
// element wise assembly
vector<int> myElems = feMesh_->owned_elts();
for (vector<int>::iterator elemsIter = myElems.begin();
elemsIter != myElems.end();
++elemsIter)
{
int ielem = *elemsIter;
// if element is masked, skip
if (elementMask && !(*elementMask)(ielem,0)) continue;
// material id
int imat = elementMaterials(ielem);
const Material * mat = physicsModel->material(imat);
// interpolate fields
interpolate_fields(ielem,fields,_conn_,_N_,_weights_,_fieldsAtIPs_);
for (int j = 0; j < nfields; ++j) {
_fieldName_ = field_mask(j);
SPAR_MAT & M = massMatrices[_fieldName_].set_quantity();
// skip null weakEqns by material
if (! (physicsModel->null(_fieldName_,imat)) ) {
physicsModel->weak_equation(_fieldName_)->
M_integrand(_fieldsAtIPs_, mat, capacity);
DIAG_MAT rho(column(capacity,0));
elementMassMatrix = _N_.transMat(_weights_*rho*_N_);
// assemble
for (int i = 0; i < nNodesPerElement_; ++i) {
int inode = _conn_(i);
for (int j = 0; j < nNodesPerElement_; ++j) {
int jnode = _conn_(j);
M.add(inode, jnode, elementMassMatrix(i,j));
}
}
}
}
}
for (int j = 0; j < nfields; ++j) {
_fieldName_ = field_mask(j);
SPAR_MAT & M = massMatrices[_fieldName_].set_quantity();
#ifdef ISOLATE_FE
sparse_allsum(communicator_,M);
#else
LammpsInterface::instance()->sparse_allsum(M);
#endif
}
// fix zero diagonal entries due to null material elements
for (int j = 0; j < nfields; ++j) {
_fieldName_ = field_mask(j);
SPAR_MAT & M = massMatrices[_fieldName_].set_quantity();
for (int inode = 0; inode < nNodesUnique_; ++inode) {
if (! M.has_entry(inode,inode)) {
M.set(inode,inode,1.0);
}
}
M.compress();
}
}
// -------------------------------------------------------------
// compute dimensionless consistent mass using native quadrature
// -------------------------------------------------------------
void FE_Engine::compute_mass_matrix(SPAR_MAT &massMatrix) const
{
// assemble nnodes X nnodes matrix
massMatrix.reset(nNodesUnique_,nNodesUnique_);// zero since partial fill
DENS_MAT elementMassMatrix(nNodesPerElement_,nNodesPerElement_);
vector<int> myElems = feMesh_->owned_elts();
for (vector<int>::iterator elemsIter = myElems.begin();
elemsIter != myElems.end();
++elemsIter)
{
int ielem = *elemsIter;
// evaluate shape functions
feMesh_->shape_function(ielem, _N_, _weights_);
// perform quadrature
elementMassMatrix = _N_.transMat(_weights_*_N_);
// get connectivity
_conn_ = feMesh_->element_connectivity_unique(ielem);
for (int i = 0; i < nNodesPerElement_; ++i) {
int inode = _conn_(i);
for (int j = 0; j < nNodesPerElement_; ++j) {
int jnode = _conn_(j);
massMatrix.add(inode, jnode, elementMassMatrix(i,j));
}
}
}
// Assemble partial results
#ifdef ISOLATE_FE
sparse_allsum(communicator_,massMatrix);
#else
LammpsInterface::instance()->sparse_allsum(massMatrix);
#endif
}
// -------------------------------------------------------------
// compute dimensionless consistent mass using given quadrature
// -------------------------------------------------------------
void FE_Engine::compute_mass_matrix(const DIAG_MAT &weights,
const SPAR_MAT &N,
SPAR_MAT &massMatrix) const
{
int nn = N.nCols();
int nips = N.nRows();
DENS_MAT tmp_mass_matrix_local(nn,nn), tmp_mass_matrix(nn,nn);
if (nips>0) { tmp_mass_matrix_local = N.transMat(weights*N); }
// share information between processors
int count = nn*nn;
allsum(communicator_,
tmp_mass_matrix_local.ptr(),
tmp_mass_matrix.ptr(), count);
// create sparse from dense
massMatrix.reset(tmp_mass_matrix,tol_sparse);
}
// -------------------------------------------------------------
// compute dimensionless lumped mass using native quadrature
// -------------------------------------------------------------
void FE_Engine::compute_lumped_mass_matrix(DIAG_MAT & M) const
{
M.reset(nNodesUnique_,nNodesUnique_); // zero since partial fill
// assemble lumped diagonal mass
DENS_VEC Nvec(nNodesPerElement_);
vector<int> myElems = feMesh_->owned_elts();
for (vector<int>::iterator elemsIter = myElems.begin();
elemsIter != myElems.end();
++elemsIter)
{
int ielem = *elemsIter;
// evaluate shape functions
feMesh_->shape_function(ielem, _N_, _weights_);
CLON_VEC w(_weights_);
// get connectivity
_conn_ = feMesh_->element_connectivity_unique(ielem);
Nvec = w*_N_;
for (int i = 0; i < nNodesPerElement_; ++i) {
int inode = _conn_(i);
M(inode,inode) += Nvec(i);
}
}
// Assemble partial results
allsum(communicator_,MPI_IN_PLACE, M.ptr(), M.size());
}
// -------------------------------------------------------------
// compute physical lumped mass using native quadrature
// -------------------------------------------------------------
void FE_Engine::compute_lumped_mass_matrix(
const Array<FieldName>& field_mask,
const FIELDS & fields,
const PhysicsModel * physicsModel,
const Array<int> &elementMaterials,
MASS_MATS & massMatrices, // mass matrix
const DenseMatrix<bool> *elementMask) const
{
int nfields = field_mask.size();
// zero initialize for assembly
for (int j = 0; j < nfields; ++j) {
DIAG_MAT & M = massMatrices[field_mask(j)].set_quantity();
M.reset(nNodesUnique_,nNodesUnique_);
}
// assemble diagonal matrix
vector<int> myElems = feMesh_->owned_elts();
for (vector<int>::iterator elemsIter = myElems.begin();
elemsIter != myElems.end();
++elemsIter)
{
int ielem = *elemsIter;
// if element is masked, skip it
if (elementMask && !(*elementMask)(ielem,0)) continue;
// material id
int imat = elementMaterials(ielem);
const Material * mat = physicsModel->material(imat);
interpolate_fields(ielem,fields,_conn_,_N_,_weights_,_fieldsAtIPs_);
// compute densities, integrate & assemble
DENS_MAT capacity;
for (int j = 0; j < nfields; ++j) {
_fieldName_ = field_mask(j);
if (! physicsModel->null(_fieldName_,imat)) {
physicsModel->weak_equation(_fieldName_)->
M_integrand(_fieldsAtIPs_, mat, capacity);
_Nmat_ = _N_.transMat(_weights_*capacity);
DIAG_MAT & M = massMatrices[_fieldName_].set_quantity();
for (int i = 0; i < nNodesPerElement_; ++i) {
int inode = _conn_(i);
M(inode,inode) += _Nmat_(i,0);
}
}
}
}
// Assemble partial results
for (int j = 0; j < nfields; ++j) {
_fieldName_ = field_mask(j);
DIAG_MAT & M = massMatrices[_fieldName_].set_quantity();
allsum(communicator_,MPI_IN_PLACE, M.ptr(), M.size());
}
// fix zero diagonal entries due to null material elements
for (int inode = 0; inode < nNodesUnique_; ++inode) {
for (int j = 0; j < nfields; ++j) {
_fieldName_ = field_mask(j);
DIAG_MAT & M = massMatrices[_fieldName_].set_quantity();
if (M(inode,inode) == 0.0) {
M(inode,inode) = 1.0;
}
}
}
}
// -------------------------------------------------------------
// compute physical lumped mass using given quadrature
// -------------------------------------------------------------
void FE_Engine::compute_lumped_mass_matrix(
const Array<FieldName> &field_mask,
const FIELDS &fields,
const PhysicsModel * physicsModel,
const Array<set<int> > & pointMaterialGroups,
const DIAG_MAT &weights,
const SPAR_MAT &N,
MASS_MATS &M) const // mass matrices
{
int nips = weights.nCols();
int nfields = field_mask.size();
// initialize
map<FieldName, DIAG_MAT> M_local;
for (int j = 0; j < nfields; ++j) {
_fieldName_ = field_mask(j);
M_local[_fieldName_].reset(nNodesUnique_,nNodesUnique_);
M[_fieldName_].reset(nNodesUnique_,nNodesUnique_);
}
if (nips>0) {
// compute fields at given ips
// compute all fields since we don't the capacities dependencies
FIELD_MATS fieldsAtIPs;
for (_fieldItr_ = fields.begin(); _fieldItr_ != fields.end(); _fieldItr_++) {
_fieldName_ = _fieldItr_->first;
const DENS_MAT & field = (_fieldItr_->second).quantity();
int dof = field.nCols();
fieldsAtIPs[_fieldName_].reset(nips,dof);
fieldsAtIPs[_fieldName_] = N*field;
}
// treat single material point sets specially
int nMatls = pointMaterialGroups.size();
int atomMatls = 0;
for (int imat = 0; imat < nMatls; imat++) {
const set<int> & indices = pointMaterialGroups(imat);
if ( indices.size() > 0) atomMatls++;
}
if (atomMatls == 0) {
throw ATC_Error("no materials in atom region - atoms may have migrated to FE-only region");
}
bool singleMaterial = ( atomMatls == 1 );
if (! singleMaterial ) {
stringstream ss; ss << " WARNING: multiple materials in atomic region";
print_msg(communicator_,ss.str());
}
// setup data structures
FIELD_MATS capacities;
// evaluate physics model & compute capacity|density for requested fields
if (singleMaterial) {
int imat = 0;
const Material * mat = physicsModel->material(imat);
for (int j = 0; j < nfields; ++j) {
_fieldName_ = field_mask(j);
// mass matrix needs to be invertible so null matls have cap=1
if (physicsModel->null(_fieldName_,imat)) {
throw ATC_Error("null material not supported for atomic region (single material)");
const FIELD & f = (fields.find(_fieldName_))->second;
capacities[_fieldName_].reset(f.nRows(),f.nCols());
capacities[_fieldName_] = 1.;
}
else {
physicsModel->weak_equation(_fieldName_)->
M_integrand(fieldsAtIPs, mat, capacities[_fieldName_]);
}
}
}
else {
FIELD_MATS groupCapacities, fieldsGroup;
for (int j = 0; j < nfields; ++j) {
_fieldName_ = field_mask(j);
capacities[_fieldName_].reset(nips,1);
}
for ( int imat = 0; imat < pointMaterialGroups.size(); imat++) {
const Material * mat = physicsModel->material(imat);
const set<int> & indices = pointMaterialGroups(imat);
int npts = indices.size();
if (npts > 0) {
for (int j = 0; j < nfields; ++j) {
_fieldName_ = field_mask(j);
groupCapacities[_fieldName_].reset(npts,1);
const FIELD & f = (fields.find(_fieldName_))->second;
int ndof = f.nCols();
fieldsGroup[_fieldName_].reset(npts,ndof);
int i = 0;
for (set<int>::const_iterator iter=indices.begin();
iter != indices.end(); iter++, i++ ) {
for (int dof = 0; dof < ndof; ++dof) {
fieldsGroup[_fieldName_](i,dof)
= fieldsAtIPs[_fieldName_](*iter,dof);
}
}
}
for (int j = 0; j < nfields; ++j) {
_fieldName_ = field_mask(j);
if (physicsModel->null(_fieldName_,imat)) {
throw ATC_Error("null material not supported for atomic region (multiple materials)");
const FIELD & f = (fields.find(_fieldName_))->second;
groupCapacities[_fieldName_].reset(f.nRows(),f.nCols());
groupCapacities[_fieldName_] = 1.;
}
else {
physicsModel->weak_equation(_fieldName_)->
M_integrand(fieldsGroup, mat, groupCapacities[_fieldName_]);
}
int i = 0;
for (set<int>::const_iterator iter=indices.begin();
iter != indices.end(); iter++, i++ ) {
capacities[_fieldName_](*iter,0)
= groupCapacities[_fieldName_](i,0);
}
}
}
}
}
// integrate & assemble
for (int j = 0; j < nfields; ++j) {
_fieldName_ = field_mask(j);
M_local[_fieldName_].reset( // assume all columns same
column(N.transMat(weights*capacities[_fieldName_]),0) );
}
}
// Share information between processors
for (int j = 0; j < nfields; ++j) {
_fieldName_ = field_mask(j);
DIAG_MAT & myMassMat(M[_fieldName_].set_quantity());
int count = M_local[_fieldName_].size();
allsum(communicator_, M_local[_fieldName_].ptr(), myMassMat.ptr(), count);
}
}
//-----------------------------------------------------------------
// compute assembled average gradient evaluated at the nodes
//-----------------------------------------------------------------
void FE_Engine::compute_gradient_matrix(SPAR_MAT_VEC & grad_matrix) const
{
// zero
DENS_MAT_VEC tmp_grad_matrix(nSD_);
for (int i = 0; i < nSD_; i++) {
tmp_grad_matrix[i].reset(nNodesUnique_,nNodesUnique_);
}
// element wise assembly
Array<int> count(nNodesUnique_); count = 0;
set_quadrature(NODAL);
vector<int> myElems = feMesh_->owned_elts();
for (vector<int>::iterator elemsIter = myElems.begin();
elemsIter != myElems.end();
++elemsIter)
{
int ielem = *elemsIter;
// evaluate shape functions
feMesh_->shape_function(ielem, _N_, _dN_, _weights_);
// get connectivity
_conn_ = feMesh_->element_connectivity_unique(ielem);
for (int j = 0; j < nIPsPerElement_; ++j) {
int jnode = _conn_(j);
count(jnode) += 1;
for (int i = 0; i < nNodesPerElement_; ++i) {
int inode = _conn_(i);
for (int k = 0; k < nSD_; ++k) {
tmp_grad_matrix[k](jnode,inode) += _dN_[k](j,i);
}
}
}
}
// Assemble partial results
for (int k = 0; k < nSD_; ++k) {
allsum(communicator_,MPI_IN_PLACE, tmp_grad_matrix[k].ptr(), tmp_grad_matrix[k].size());
}
int_allsum(communicator_,MPI_IN_PLACE, count.ptr(), count.size());
set_quadrature(quadrature_); //reset to default
for (int inode = 0; inode < nNodesUnique_; ++inode) {
for (int jnode = 0; jnode < nNodesUnique_; ++jnode) {
for (int k = 0; k < nSD_; ++k) {
tmp_grad_matrix[k](jnode,inode) /= count(jnode);
}
}
}
// compact dense matrices
for (int k = 0; k < nSD_; ++k) {
grad_matrix[k]->reset(tmp_grad_matrix[k],tol_sparse);
}
}
// -------------------------------------------------------------
// compute energy per node using native quadrature
// -------------------------------------------------------------
void FE_Engine::compute_energy(const Array<FieldName> &mask,
const FIELDS &fields,
const PhysicsModel * physicsModel,
const Array<int> & elementMaterials,
FIELD_MATS &energies,
const DenseMatrix<bool> *elementMask,
const IntegrationDomainType domain) const
{
// Zero out all fields
for (int n = 0; n < mask.size(); n++) {
_fieldName_ = mask(n);
_fieldItr_ = fields.find(_fieldName_);
const DENS_MAT & field = (_fieldItr_->second).quantity();
energies[_fieldName_].reset(field.nRows(), 1);
}
DENS_MAT elementEnergy(nNodesPerElement_,1); // workspace
vector<int> myElems = feMesh_->owned_elts();
for (vector<int>::iterator elemsIter = myElems.begin();
elemsIter != myElems.end();
++elemsIter)
{
int ielem = *elemsIter;
// if element is masked, skip it
if (domain != FULL_DOMAIN && elementMask && !(*elementMask)(ielem,0)) continue;
// material id
int imat = elementMaterials(ielem);
const Material * mat = physicsModel->material(imat);
interpolate_fields(ielem,fields,_conn_,_N_,_dN_,_weights_,
_fieldsAtIPs_,_gradFieldsAtIPs_);
// assemble
for (int n = 0; n < mask.size(); n++) {
_fieldName_ = mask(n);
if (physicsModel->null(_fieldName_,imat)) continue;
if( ! (physicsModel->weak_equation(_fieldName_)-> has_E_integrand())) continue;
physicsModel->weak_equation(_fieldName_)->
E_integrand(_fieldsAtIPs_, _gradFieldsAtIPs_, mat, elementEnergy);
_fieldItr_ = fields.find(_fieldName_);
_Nmat_ = _N_.transMat(_weights_*elementEnergy);
DENS_MAT & energy = energies[_fieldName_];
for (int i = 0; i < nNodesPerElement_; ++i) {
int inode = _conn_(i);
energy(inode,0) += _Nmat_(i,0);
}
}
}
for (int n = 0; n < mask.size(); n++) {
_fieldName_ = mask(n);
DENS_MAT& myEnergy(energies[_fieldName_]);
allsum(communicator_,MPI_IN_PLACE, myEnergy.ptr(), myEnergy.size());
}
}
// -------------------------------------------------------------
// compute rhs using native quadrature
// -------------------------------------------------------------
void FE_Engine::compute_rhs_vector(const Array2D<bool> &rhsMask,
const FIELDS &fields,
const PhysicsModel * physicsModel,
const Array<int> & elementMaterials,
FIELDS &rhs,
const DenseMatrix<bool> *elementMask) const
{
vector<FieldName> usedFields;
// size and zero output
for (int n = 0; n < rhsMask.nRows(); n++) {
_fieldName_ = FieldName(n);
if (rhsMask(_fieldName_, FLUX) || rhsMask(_fieldName_, SOURCE)) {
_fieldItr_ = fields.find(_fieldName_);
const DENS_MAT & field = (_fieldItr_->second).quantity();
rhs[_fieldName_].reset(field.nRows(), field.nCols());
// Save field names for easy lookup later.
usedFields.push_back(_fieldName_);
}
}
// Iterate over elements partitioned onto current processor.
vector<int> myElems = feMesh_->owned_elts();
for (vector<int>::iterator elemsIter = myElems.begin();
elemsIter != myElems.end();
++elemsIter)
{
int ielem = *elemsIter;
// Skip masked elements
if (elementMask && !(*elementMask)(ielem,0)) continue;
int imat = elementMaterials(ielem); // material id
const Material * mat = physicsModel->material(imat);
// interpolate field values to integration points
interpolate_fields(ielem,fields,_conn_,_N_,_dN_,_weights_,
_fieldsAtIPs_,_gradFieldsAtIPs_);
// rescale by _weights_, a diagonal matrix
_Nw_ = _weights_*_N_;
for (int j = 0; j < nSD_; ++j) _dNw_[j] = _weights_*_dN_[j];
// evaluate physics model and assemble
// _Nfluxes is a set of fields
for (int n = 0; n < rhsMask.nRows(); n++) {
_fieldName_ = FieldName(n);
if (!rhsMask(_fieldName_,SOURCE)) continue;
if (physicsModel->null(_fieldName_,imat)) continue;
_fieldItr_ = fields.find(_fieldName_);
const DENS_MAT & field = (_fieldItr_->second).quantity();
DENS_MAT & myRhs(rhs[_fieldName_].set_quantity());
int dof = field.nCols();
bool has = physicsModel->weak_equation(_fieldName_)->
N_integrand(_fieldsAtIPs_,_gradFieldsAtIPs_, mat, _Nfluxes_);
if (!has) continue;
_Nmat_ = _Nw_.transMat(_Nfluxes_);
for (int i = 0; i < nNodesPerElement_; ++i) {
int inode = _conn_(i);
for (int k = 0; k < dof; ++k) {
myRhs(inode,k) += _Nmat_(i,k);
}
}
}
// _Bfluxes_ is a set of field gradients
for (int n = 0; n < rhsMask.nRows(); n++) {
_fieldName_ = FieldName(n);
if (!rhsMask(_fieldName_,FLUX)) continue;
if (physicsModel->null(_fieldName_,imat)) continue;
_fieldItr_ = fields.find(_fieldName_);
const DENS_MAT & field = (_fieldItr_->second).quantity();
DENS_MAT & myRhs(rhs[_fieldName_].set_quantity());
int dof = field.nCols();
physicsModel->weak_equation(_fieldName_)->
B_integrand(_fieldsAtIPs_, _gradFieldsAtIPs_, mat, _Bfluxes_);
for (int j = 0; j < nSD_; j++) {
_Nmat_ = _dNw_[j].transMat(_Bfluxes_[j]);
for (int i = 0; i < nNodesPerElement_; ++i) {
int inode = _conn_(i);
for (int k = 0; k < dof; ++k) {
myRhs(inode,k) += _Nmat_(i,k);
}
}
}
}
}
vector<FieldName>::iterator _fieldIter_;
for (_fieldIter_ = usedFields.begin(); _fieldIter_ != usedFields.end();
++_fieldIter_) {
// myRhs is where we put the global result for this field.
_fieldName_ = *_fieldIter_;
DENS_MAT & myRhs(rhs[_fieldName_].set_quantity());
// Sum matrices in-place across all processors into myRhs.
allsum(communicator_,MPI_IN_PLACE, myRhs.ptr(), myRhs.size());
}
}
// -------------------------------------------------------------
// compute rhs using given quadrature
// -------------------------------------------------------------
void FE_Engine::compute_rhs_vector(const Array2D<bool> &rhsMask,
const FIELDS &fields,
const PhysicsModel *physicsModel,
const Array<set<int> >&pointMaterialGroups,
const DIAG_MAT &weights,
const SPAR_MAT &N,
const SPAR_MAT_VEC &dN,
FIELDS &rhs) const
{
FIELD_MATS rhs_local;
for (int n = 0; n < rhsMask.nRows(); n++) {
_fieldName_ = FieldName(n);
if (rhsMask(_fieldName_,FLUX) || rhsMask(_fieldName_,SOURCE)) {
_fieldItr_ = fields.find(_fieldName_);
const DENS_MAT & field = (_fieldItr_->second).quantity();
int nrows = field.nRows();
int ncols = field.nCols();
rhs [_fieldName_].reset(nrows,ncols);
rhs_local[_fieldName_].reset(nrows,ncols);
}
}
int nips = weights.nCols();
if (nips>0) {
// compute fields and gradients of fields at given ips
GRAD_FIELD_MATS gradFieldsAtIPs;
FIELD_MATS fieldsAtIPs;
for (_fieldItr_ = fields.begin(); _fieldItr_ != fields.end(); _fieldItr_++) {
_fieldName_ = _fieldItr_->first;
const DENS_MAT & field = (_fieldItr_->second).quantity();
int dof = field.nCols();
gradFieldsAtIPs[_fieldName_].assign(nSD_,DENS_MAT(nips,dof));
fieldsAtIPs[_fieldName_] = N*field;
for (int j = 0; j < nSD_; ++j) {
gradFieldsAtIPs[_fieldName_][j] = (*dN[j])*field;
}
}
// setup data structures
FIELD_MATS Nfluxes;
GRAD_FIELD_MATS Bfluxes;
for (int n = 0; n < rhsMask.nRows(); n++) {
_fieldName_ = FieldName(n);
int index = (int) _fieldName_;
if ( rhsMask(index,FLUX) ) {
Bfluxes[_fieldName_].assign(nSD_, DENS_MAT());
}
}
// treat single material point sets specially
int nMatls = pointMaterialGroups.size();
int atomMatls = 0;
for (int imat = 0; imat < nMatls; imat++) {
const set<int> & indices = pointMaterialGroups(imat);
if ( indices.size() > 0) atomMatls++;
}
bool singleMaterial = ( atomMatls == 1 );
// evaluate physics model
if (singleMaterial)
{
int imat = 0;
const Material * mat = physicsModel->material(imat);
for (int n = 0; n < rhsMask.nRows(); n++) {
_fieldName_ = FieldName(n);
int index = (int) _fieldName_;
if (! physicsModel->null(_fieldName_,imat)) {
if ( rhsMask(index,SOURCE) ) {
physicsModel->weak_equation(_fieldName_)->
N_integrand(fieldsAtIPs, gradFieldsAtIPs, mat, Nfluxes[_fieldName_]);
}
if ( rhsMask(index,FLUX) ) {
physicsModel->weak_equation(_fieldName_)->
B_integrand(fieldsAtIPs, gradFieldsAtIPs, mat, Bfluxes[_fieldName_]);
}
}
}
}
else
{
// 1) permanent workspace with per-row mapped clones per matl
// from caller/atc
// 2) per point calls to N/B_integrand
// 3) collect internalToAtom into materials and use mem-cont clones
// what about dof for matrices and data storage: clone _matrix_
// for each material group:
// set up storage
DENS_MAT group_Nfluxes;
DENS_MAT_VEC group_Bfluxes;
group_Bfluxes.reserve(nSD_);
FIELD_MATS fieldsGroup;
GRAD_FIELD_MATS gradFieldsGroup;
for (_fieldItr_ = fields.begin(); _fieldItr_ != fields.end(); _fieldItr_++) {
_fieldName_ = _fieldItr_->first;
const DENS_MAT & field = (_fieldItr_->second).quantity();
int ndof = field.nCols();
gradFieldsGroup[_fieldName_].assign(nSD_,DENS_MAT(nips,ndof));
Nfluxes[_fieldName_].reset(nips,ndof);
Bfluxes[_fieldName_].assign(nSD_,DENS_MAT(nips,ndof));
//}
}
// copy fields
for ( int imat = 0; imat < pointMaterialGroups.size(); imat++)
{
const set<int> & indices = pointMaterialGroups(imat);
const Material * mat = physicsModel->material(0);
int npts = indices.size();
int i = 0;
for (set<int>::const_iterator iter=indices.begin();
iter != indices.end(); iter++, i++ ) {
for (_fieldItr_ = fields.begin(); _fieldItr_ != fields.end(); _fieldItr_++) {
_fieldName_ = _fieldItr_->first;
const DENS_MAT & field = (_fieldItr_->second).quantity();
int ndof = field.nCols();
fieldsGroup[_fieldName_].reset(npts,ndof);
for (int j = 0; j < nSD_; ++j) {
(gradFieldsGroup[_fieldName_][j]).reset(npts,ndof);
}
for (int dof = 0; dof < ndof; ++dof) {
fieldsGroup[_fieldName_](i,dof)
= fieldsAtIPs[_fieldName_](*iter,dof);
for (int j = 0; j < nSD_; ++j) {
gradFieldsGroup[_fieldName_][j](i,dof)
= gradFieldsAtIPs[_fieldName_][j](*iter,dof);
}
}
}
}
// calculate forces, & assemble
for (int n = 0; n < rhsMask.nRows(); n++) {
_fieldName_ = FieldName(n);
_fieldItr_ = fields.find(_fieldName_);
int index = (int) _fieldName_;
if (physicsModel->null(_fieldName_,imat)) continue;
if ( rhsMask(index,SOURCE) ) {
const DENS_MAT & field = (_fieldItr_->second).quantity();
int ndof = field.nCols();
bool has = physicsModel->weak_equation(_fieldName_)->
N_integrand(fieldsGroup, gradFieldsGroup, mat, group_Nfluxes);
if (! has) throw ATC_Error("atomic source can not be null currently");
int i = 0;
for (set<int>::const_iterator iter=indices.begin();
iter != indices.end(); iter++, i++ ) {
for (int dof = 0; dof < ndof; ++dof) {
Nfluxes[_fieldName_](*iter,dof) += group_Nfluxes(i,dof);
}
}
}
}
// calculate forces, & assemble
for (int n = 0; n < rhsMask.nRows(); n++) {
_fieldName_ = FieldName(n);
if (physicsModel->null(_fieldName_,imat)) continue;
int index = (int) _fieldName_;
if ( rhsMask(index,FLUX) ) {
_fieldItr_ = fields.find(_fieldName_);
const DENS_MAT & field = (_fieldItr_->second).quantity();
int ndof = field.nCols();
physicsModel->weak_equation(_fieldName_)->
B_integrand(fieldsGroup, gradFieldsGroup, mat, group_Bfluxes);
int i = 0;
for (set<int>::const_iterator iter=indices.begin();
iter != indices.end(); iter++, i++ ) {
for (int dof = 0; dof < ndof; ++dof) {
for (int j = 0; j < nSD_; ++j) {
Bfluxes[_fieldName_][j](*iter,dof) += group_Bfluxes[j](i,dof);
}
}
}
}
}
}
} // endif multiple materials
// assemble : N/Bfluxes -> rhs_local
for (int n = 0; n < rhsMask.nRows(); n++) {
_fieldName_ = FieldName(n);
int index = (int) _fieldName_;
if ( rhsMask(index,SOURCE) && Nfluxes[_fieldName_].nCols() > 0 ) {
rhs_local[_fieldName_]
+= N.transMat(weights*Nfluxes[_fieldName_]);
}
if ( rhsMask(index,FLUX) && (Bfluxes[_fieldName_][0]).nCols() > 0 ) {
for (int j = 0; j < nSD_; ++j) {
rhs_local[_fieldName_]
+= dN[j]->transMat(weights*Bfluxes[_fieldName_][j]);
}
}
}
} // end nips > 0
// Share information between processors: rhs_local -> rhs
for (int n = 0; n < rhsMask.nRows(); n++) {
_fieldName_ = FieldName(n);
if (rhsMask(_fieldName_,FLUX) || rhsMask(_fieldName_,SOURCE)) {
DENS_MAT & myRhs(rhs[_fieldName_].set_quantity());
int count = rhs_local[_fieldName_].size();
allsum(communicator_, rhs_local[_fieldName_].ptr(), myRhs.ptr(), count);
}
}
}
// -------------------------------------------------------------
// compute sources using given quadrature
// -------------------------------------------------------------
void FE_Engine::compute_source(const Array2D<bool> &rhsMask,
const FIELDS &fields,
const PhysicsModel *physicsModel,
const Array<set<int> >&pointMaterialGroups,
const DIAG_MAT &weights,
const SPAR_MAT &N,
const SPAR_MAT_VEC &dN,
FIELD_MATS &sources) const
{
int nips = weights.nCols();
if (nips>0) {
FIELD_MATS Nfluxes;
// compute fields and gradients of fields at given ips
GRAD_FIELD_MATS gradFieldsAtIPs;
FIELD_MATS fieldsAtIPs;
for (_fieldItr_ = fields.begin(); _fieldItr_ != fields.end(); _fieldItr_++) {
_fieldName_ = _fieldItr_->first;
const DENS_MAT & field = (_fieldItr_->second).quantity();
int dof = field.nCols();
gradFieldsAtIPs[_fieldName_].assign(nSD_,DENS_MAT(nips,dof));
fieldsAtIPs[_fieldName_] = N*field;
for (int j = 0; j < nSD_; ++j) {
gradFieldsAtIPs[_fieldName_][j] = (*dN[j])*field;
}
}
// treat single material point sets specially
int nMatls = pointMaterialGroups.size();
int atomMatls = 0;
for (int imat = 0; imat < nMatls; imat++) {
const set<int> & indices = pointMaterialGroups(imat);
if ( indices.size() > 0) atomMatls++;
}
bool singleMaterial = ( atomMatls == 1 );
// evaluate physics model
if (singleMaterial)
{
int imat = 0;
const Material * mat = physicsModel->material(imat);
for (int n = 0; n < rhsMask.nRows(); n++) {
_fieldName_ = FieldName(n);
int index = (int) _fieldName_;
if (! physicsModel->null(_fieldName_,imat)) {
if ( rhsMask(index,SOURCE) ) {
bool has = physicsModel->weak_equation(_fieldName_)->
N_integrand(fieldsAtIPs, gradFieldsAtIPs, mat, Nfluxes[_fieldName_]);
if (! has) throw ATC_Error("atomic source can not be null currently");
}
}
}
}
else
{
throw ATC_Error("compute source does not handle multiple materials currently");
}
// assemble
for (int n = 0; n < rhsMask.nRows(); n++) {
_fieldName_ = FieldName(n);
int index = (int) _fieldName_;
if ( rhsMask(index,SOURCE) ) {
sources[_fieldName_] =weights*Nfluxes[_fieldName_];
}
}
}
// no need to share information between processors
}
// -------------------------------------------------------------
// compute flux for post processing
// -------------------------------------------------------------
void FE_Engine::compute_flux(const Array2D<bool> &rhsMask,
const FIELDS &fields,
const PhysicsModel * physicsModel,
const Array<int> & elementMaterials,
GRAD_FIELD_MATS &fluxes,
const DenseMatrix<bool> *elementMask) const
{
for (int n = 0; n < rhsMask.nRows(); n++) {
_fieldName_ = FieldName(n);
if (rhsMask(_fieldName_,FLUX)) {
_fieldItr_ = fields.find(_fieldName_);
const DENS_MAT & field = (_fieldItr_->second).quantity();
fluxes[_fieldName_].assign(nSD_, DENS_MAT(field.nRows(),field.nCols()));
}
}
// element wise assembly
vector<int> myElems = feMesh_->owned_elts();
for (vector<int>::iterator elemsIter = myElems.begin();
elemsIter != myElems.end();
++elemsIter)
{
int ielem = *elemsIter;
// if element is masked, skip it
if (elementMask && !(*elementMask)(ielem,0)) continue;
// material id
int imat = elementMaterials(ielem);
const Material * mat = physicsModel->material(imat);
interpolate_fields(ielem,fields,_conn_,_N_,_dN_,_weights_,
_fieldsAtIPs_,_gradFieldsAtIPs_);
_Nw_ = _weights_*_N_;
// evaluate Physics model & assemble
for (int n = 0; n < rhsMask.nRows(); n++) {
_fieldName_ = FieldName(n);
if (!rhsMask(_fieldName_,FLUX)) continue;
if (physicsModel->null(_fieldName_,imat)) continue;
_fieldItr_ = fields.find(_fieldName_);
const DENS_MAT & field = (_fieldItr_->second).quantity();
int dof = field.nCols();
physicsModel->weak_equation(_fieldName_)->
B_integrand(_fieldsAtIPs_, _gradFieldsAtIPs_, mat, _Bfluxes_);
for (int j = 0; j < nSD_; j++) {
_Nmat_ = _Nw_.transMat(_Bfluxes_[j]);
for (int i = 0; i < nNodesPerElement_; ++i) {
int inode = _conn_(i);
for (int k = 0; k < dof; ++k) {
fluxes[_fieldName_][j](inode,k) += _Nmat_(i,k);
}
}
}
}
}
// Assemble partial results
for (int n = 0; n < rhsMask.nRows(); n++) {
_fieldName_ = FieldName(n);
if (!rhsMask(_fieldName_,FLUX)) continue;
for (int j = 0; j < nSD_; j++) {
allsum(communicator_,MPI_IN_PLACE, fluxes[_fieldName_][j].ptr(), fluxes[_fieldName_][j].size());
}
}
}
//-----------------------------------------------------------------
// boundary flux using native quadrature
//-----------------------------------------------------------------
void FE_Engine::compute_boundary_flux(const Array2D<bool> & rhsMask,
const FIELDS & fields,
const PhysicsModel * physicsModel,
const Array<int> & elementMaterials,
const set< pair <int,int> > & faceSet,
FIELDS & rhs) const
{
for (int n = 0; n < rhsMask.nRows(); n++) {
_fieldName_ = FieldName(n);
if (rhsMask(_fieldName_,FLUX)) {
_fieldItr_ = fields.find(_fieldName_);
const DENS_MAT & field = (_fieldItr_->second).quantity();
rhs[_fieldName_].reset(field.nRows(),field.nCols());
}
}
FIELD_MATS localElementFields;
GRAD_FIELD_MATS gradFieldsAtIPs;
FIELD_MATS fieldsAtIPs;
for (_fieldItr_ = fields.begin(); _fieldItr_ != fields.end(); _fieldItr_++) {
_fieldName_ = _fieldItr_->first;
const FIELD & field = _fieldItr_->second;
int dof = field.nCols();
gradFieldsAtIPs[_fieldName_].reserve(nSD_);
for (int i = 0; i < nSD_; ++i) {
gradFieldsAtIPs[_fieldName_].push_back(DENS_MAT(nIPsPerFace_,dof));
}
fieldsAtIPs[_fieldName_].reset(nIPsPerFace_,dof);
localElementFields[_fieldName_].reset(nNodesPerElement_,dof);
}
// element wise assembly
set< PAIR >::iterator iter;
for (iter = faceSet.begin(); iter != faceSet.end(); iter++) {
// get connectivity
int ielem = iter->first;
// if this is not our element, do not do calculations
if (!feMesh_->is_owned_elt(ielem)) continue;
int imat = elementMaterials(ielem);
const Material * mat = physicsModel->material(imat);
_conn_ = feMesh_->element_connectivity_unique(ielem);
// evaluate shape functions at ips
feMesh_->face_shape_function(*iter, _fN_, _fdN_, _nN_, _fweights_);
// interpolate fields and gradients of fields ips of this element
for (_fieldItr_ = fields.begin(); _fieldItr_ != fields.end(); _fieldItr_++) {
_fieldName_ = _fieldItr_->first;
const DENS_MAT & field = (_fieldItr_->second).quantity();
int dof = field.nCols();
for (int i = 0; i < nNodesPerElement_; i++) {
for (int j = 0; j < dof; j++) {
localElementFields[_fieldName_](i,j) = field(_conn_(i),j);
}
}
// ips X dof = ips X nodes * nodes X dof
fieldsAtIPs[_fieldName_] = _fN_*localElementFields[_fieldName_];
for (int j = 0; j < nSD_; ++j) {
gradFieldsAtIPs[_fieldName_][j] = _fdN_[j]*localElementFields[_fieldName_];
}
}
// Evaluate- physics model
// do nothing for N_integrand
// nN : precomputed and held by ATC_Transfer
// assemble
for (int n = 0; n < rhsMask.nRows(); n++) {
_fieldName_ = FieldName(n);
int index = (int) _fieldName_;
if ( rhsMask(index,FLUX) ) {
_fieldItr_ = fields.find(_fieldName_);
const DENS_MAT & field = (_fieldItr_->second).quantity();
DENS_MAT & myRhs(rhs[_fieldName_].set_quantity());
physicsModel->weak_equation(_fieldName_)->
B_integrand(fieldsAtIPs, gradFieldsAtIPs, mat, _Bfluxes_);
int dof = field.nCols();
for (int j = 0; j < nSD_; j++) {
_Nmat_ = _nN_[j].transMat(_fweights_*_Bfluxes_[j]);
for (int i = 0; i < nNodesPerElement_; ++i) {
int inode = _conn_(i);
for (int k = 0; k < dof; ++k) {
myRhs(inode,k) += _Nmat_(i,k);
}
}
}
}
}
}
for (int n = 0; n < rhsMask.nRows(); n++) {
_fieldName_ = FieldName(n);
int index = (int) _fieldName_;
if ( rhsMask(index,FLUX) ) {
DENS_MAT & myRhs(rhs[_fieldName_].set_quantity());
allsum(communicator_,MPI_IN_PLACE, myRhs.ptr(), myRhs.size());
}
}
}
// -------------------------------------------------------------
// compute boundary flux using given quadrature and interpolation
// -------------------------------------------------------------
void FE_Engine::compute_boundary_flux(const Array2D<bool> & rhsMask,
const FIELDS & fields,
const PhysicsModel * physicsModel,
const Array< int > & elementMaterials,
const Array< set<int> > & pointMaterialGroups,
const DIAG_MAT & _weights_,
const SPAR_MAT & N,
const SPAR_MAT_VEC & dN,
const DIAG_MAT & flux_mask,
FIELDS & flux,
const DenseMatrix<bool> * elementMask,
const set<int> * nodeSet) const
{
FIELDS rhs, rhsSubset;
for (int n = 0; n < rhsMask.nRows(); n++) {
_fieldName_ = FieldName(n);
if (rhsMask(_fieldName_,FLUX)) {
_fieldItr_ = fields.find(_fieldName_);
const DENS_MAT & field = (_fieldItr_->second).quantity();
int nrows = field.nRows();
int ncols = field.nCols();
rhs [_fieldName_].reset(nrows,ncols);
rhsSubset[_fieldName_].reset(nrows,ncols);
}
}
// R_I = - int_Omega Delta _N_I .q dV
compute_rhs_vector(rhsMask, fields, physicsModel, elementMaterials, rhs, elementMask);
// R_I^md = - int_Omega^md Delta _N_I .q dV
compute_rhs_vector(rhsMask, fields, physicsModel, pointMaterialGroups,
_weights_, N, dN, rhsSubset);
// flux_I = 1/Delta V_I V_I^md R_I + R_I^md
// where : Delta V_I = int_Omega _N_I dV
for (int n = 0; n < rhsMask.nRows(); n++) {
_fieldName_ = FieldName(n);
if ( rhsMask(_fieldName_,FLUX) ) {
if (nodeSet) {
_fieldItr_ = fields.find(_fieldName_);
const DENS_MAT & field = (_fieldItr_->second).quantity();
int nrows = field.nRows();
int ncols = field.nCols();
DENS_MAT & myFlux(flux[_fieldName_].set_quantity());
const DENS_MAT & myRhsSubset(rhsSubset[_fieldName_].quantity());
const DENS_MAT & myRhs(rhs[_fieldName_].quantity());
myFlux.reset(nrows,ncols);
set<int>::const_iterator iset;
for (iset = nodeSet->begin(); iset != nodeSet->end(); iset++) {
for (int j = 0; j < ncols; j++) {
myFlux(*iset,j) = myRhsSubset(*iset,j) - flux_mask(*iset,*iset)*myRhs(*iset,j);
}
}
}
else {
flux[_fieldName_]
= rhsSubset[_fieldName_].quantity() - flux_mask*(rhs[_fieldName_].quantity());
}
}
}
}
/** integrate a nodal field over an element set */
DENS_VEC FE_Engine::integrate(const DENS_MAT &field, const ESET & eset) const
{
int dof = field.nCols();
DENS_MAT eField(nNodesPerElement_, dof);
int nips = nIPsPerElement_;
DENS_MAT ipField(nips, dof);
DENS_VEC integral(dof); integral = 0;
for (ESET::const_iterator itr = eset.begin(); itr != eset.end(); ++itr) {
int ielem = *itr;
// if this is not our element, do not do calculations
if (!feMesh_->is_owned_elt(ielem)) continue;
feMesh_->shape_function(ielem,_N_, _weights_);
_conn_ = feMesh_->element_connectivity_unique(ielem);
for (int i = 0; i < nNodesPerElement_; i++) {
for (int j = 0; j < dof; j++) {
eField(i,j) = field(_conn_(i),j); }}
ipField = _N_*eField;
for (int i = 0; i < nips; ++i) {
for (int j = 0; j < dof; ++j) {
integral(j) += ipField(i,j)*_weights_[i];
}
}
}
// assemble partial results
allsum(communicator_,MPI_IN_PLACE, integral.ptr(), integral.size());
return integral;
}
//-----------------------------------------------------------------
// Robin boundary flux using native quadrature
// integrate \int_delV _N_I s(x,t).n dA
//-----------------------------------------------------------------
void FE_Engine::add_robin_fluxes(const Array2D<bool> &rhsMask,
const FIELDS & fields,
const double time,
const ROBIN_SURFACE_SOURCE & sourceFunctions,
FIELDS &nodalSources) const
{
// sizing working arrays
DENS_MAT xCoords(nSD_,nNodesPerElement_);
DENS_MAT faceSource;
DENS_MAT localField;
DENS_MAT xAtIPs(nSD_,nIPsPerFace_);
DENS_MAT uAtIPs(nIPsPerFace_,1);
// element wise assembly
ROBIN_SURFACE_SOURCE::const_iterator src_iter;
if (!(rank_zero(communicator_))) {
// Zero out unmasked nodal sources on all non-main processors.
// This is to avoid counting the previous nodal source values
// multiple times when we aggregate.
for (src_iter=sourceFunctions.begin(); src_iter!=sourceFunctions.end(); src_iter++)
{
_fieldName_ = src_iter->first;
if (!rhsMask((int)_fieldName_,ROBIN_SOURCE)) continue;
DENS_MAT & myNodalSources(nodalSources[_fieldName_].set_quantity());
myNodalSources.reset(myNodalSources.nRows(), myNodalSources.nCols());
}
}
for (src_iter=sourceFunctions.begin(); src_iter!=sourceFunctions.end(); src_iter++)
{
_fieldName_ = src_iter->first;
if (!rhsMask((int)_fieldName_,ROBIN_SOURCE)) continue;
typedef map<PAIR,Array<UXT_Function*> > FSET;
const FSET *fset = (const FSET *)&(src_iter->second);
FSET::const_iterator fset_iter;
for (fset_iter = fset->begin(); fset_iter != fset->end(); fset_iter++)
{
const PAIR &face = fset_iter->first;
const int elem = face.first;
// if this is not our element, do not do calculations
if (!feMesh_->is_owned_elt(elem)) continue;
const Array <UXT_Function*> &fs = fset_iter->second;
_conn_ = feMesh_->element_connectivity_unique(elem);
// evaluate location at ips
feMesh_->face_shape_function(face, _fN_, _fdN_, _nN_, _fweights_);
feMesh_->element_coordinates(elem, xCoords);
xAtIPs = xCoords*(_fN_.transpose());
// collect field
_fieldItr_ = fields.find(_fieldName_);
const DENS_MAT & field = (_fieldItr_->second).quantity();
feMesh_->element_field(elem, field, localField);
uAtIPs = _fN_*localField;
// interpolate prescribed flux at ips of this element
FSET::const_iterator face_iter = fset->find(face);
int nFieldDOF = (face_iter->second).size();
faceSource.reset(nIPsPerFace_,nFieldDOF);
for (int ip = 0; ip < nIPsPerFace_; ++ip) {
for (int idof = 0; idof<nFieldDOF; ++idof) {
UXT_Function * f = fs(idof);
if (!f) continue;
faceSource(ip,idof) = f->f(&(uAtIPs(ip,0)),
column(xAtIPs,ip).ptr(),time);
DENS_MAT coefsAtIPs(nIPsPerFace_,1);
coefsAtIPs(ip,idof) = f->dfdu(&(uAtIPs(ip,0)),
column(xAtIPs,ip).ptr(),time);
faceSource(ip,idof) -= coefsAtIPs(ip,0)*uAtIPs(ip,0);
}
}
// assemble
DENS_MAT & myNodalSources(nodalSources[_fieldName_].set_quantity());
_Nmat_ = _fN_.transMat(_fweights_*faceSource);
for (int i = 0; i < nNodesPerElement_; ++i) {
int inode = _conn_(i);
for (int idof = 0; idof < nFieldDOF; ++idof) {
myNodalSources(inode,idof) += _Nmat_(i,idof);
}
}
}
}
// assemble partial result matrices
for (src_iter=sourceFunctions.begin(); src_iter!=sourceFunctions.end(); src_iter++) {
_fieldName_ = src_iter->first;
if (!rhsMask((int) _fieldName_,ROBIN_SOURCE)) continue;
DENS_MAT & myNodalSources(nodalSources[_fieldName_].set_quantity());
allsum(communicator_,MPI_IN_PLACE, myNodalSources.ptr(), myNodalSources.size());
}
}
//-----------------------------------------------------------------
// Robin boundary flux stiffness using native quadrature
// integrate \int_delV _N_I ds/du(x,t).n dA
//-----------------------------------------------------------------
void FE_Engine::add_robin_tangent(const Array2D<bool> &rhsMask,
const FIELDS & fields,
const double time,
const ROBIN_SURFACE_SOURCE & sourceFunctions,
SPAR_MAT &tangent) const
{
// sizing working arrays
DENS_MAT xCoords(nSD_,nNodesPerElement_);
DENS_MAT coefsAtIPs;
DENS_MAT localField;
DENS_MAT xAtIPs(nSD_,nIPsPerFace_);
DENS_MAT uAtIPs(nIPsPerFace_,1);
// element wise assembly
ROBIN_SURFACE_SOURCE::const_iterator src_iter;
if (!(rank_zero(communicator_))) {
// Zero out result (tangent) matrix on all non-main processors
// to avoid multiple-counting of the values already in tangent
tangent.reset(tangent.nRows(), tangent.nCols());
}
for (src_iter=sourceFunctions.begin(); src_iter!=sourceFunctions.end(); src_iter++)
{
_fieldName_ = src_iter->first;
if (!rhsMask((int)_fieldName_,ROBIN_SOURCE)) continue;
typedef map<PAIR,Array<UXT_Function*> > FSET;
const FSET *fset = (const FSET *)&(src_iter->second);
FSET::const_iterator fset_iter;
for (fset_iter = fset->begin(); fset_iter != fset->end(); fset_iter++)
{
const PAIR &face = fset_iter->first;
const int elem = face.first;
// if this is not our element, do not do calculations
if (!feMesh_->is_owned_elt(elem)) continue;
const Array <UXT_Function*> &fs = fset_iter->second;
_conn_ = feMesh_->element_connectivity_unique(elem);
// evaluate location at ips
feMesh_->face_shape_function(face, _fN_, _fdN_, _nN_, _fweights_);
feMesh_->element_coordinates(elem, xCoords);
xAtIPs = xCoords*(_fN_.transpose());
// collect field
_fieldItr_ = fields.find(_fieldName_);
const DENS_MAT & field = (_fieldItr_->second).quantity();
feMesh_->element_field(elem, field, localField);
uAtIPs = _fN_*localField;
// interpolate prescribed flux at ips of this element
FSET::const_iterator face_iter = fset->find(face);
int nFieldDOF = (face_iter->second).size();
coefsAtIPs.reset(nIPsPerFace_,nFieldDOF);
for (int ip = 0; ip < nIPsPerFace_; ++ip) {
for (int idof = 0; idof<nFieldDOF; ++idof) {
UXT_Function * f = fs(idof);
if (!f) continue;
coefsAtIPs(ip,idof) = f->dfdu(&(uAtIPs(ip,0)),
column(xAtIPs,ip).ptr(),time);
}
}
// assemble
DIAG_MAT D(column(coefsAtIPs,0));
D *= -1;
D *= _fweights_;
_Nmat_ = _fN_.transMat(D*_fN_);
for (int i = 0; i < nNodesPerElement_; ++i) {
int inode = _conn_(i);
for (int j = 0; j < nNodesPerElement_; ++j) {
int jnode = _conn_(j);
tangent.add(inode, jnode, _Nmat_(i,j));
}
}
}
}
// assemble partial result matrices
#ifdef ISOLATE_FE
sparse_allsum(communicator_,tangent);
#else
LammpsInterface::instance()->sparse_allsum(tangent);
#endif
}
//-----------------------------------------------------------------
// prescribed boundary flux using native quadrature
// integrate \int_delV _N_I s(x,t).n dA
//-----------------------------------------------------------------
void FE_Engine::add_fluxes(const Array<bool> &fieldMask,
const double time,
const SURFACE_SOURCE & sourceFunctions,
FIELDS &nodalSources) const
{
// sizing working arrays
DENS_MAT xCoords(nSD_,nNodesPerElement_);
DENS_MAT xAtIPs(nSD_,nIPsPerFace_);
DENS_MAT faceSource;
// element wise assembly
SURFACE_SOURCE::const_iterator src_iter;
if (!(rank_zero(communicator_))) {
// Zero out unmasked nodal sources on all non-main processors.
// This is to avoid counting the previous nodal source values
// multiple times when we aggregate.
for (src_iter=sourceFunctions.begin(); src_iter!=sourceFunctions.end(); src_iter++)
{
_fieldName_ = src_iter->first;
if (!fieldMask((int)_fieldName_)) continue;
DENS_MAT & myNodalSources(nodalSources[_fieldName_].set_quantity());
myNodalSources.reset(myNodalSources.nRows(), myNodalSources.nCols());
}
}
for (src_iter=sourceFunctions.begin(); src_iter!=sourceFunctions.end(); src_iter++)
{
_fieldName_ = src_iter->first;
if (!fieldMask((int)_fieldName_)) continue;
typedef map<PAIR,Array<XT_Function*> > FSET;
const FSET *fset = (const FSET *)&(src_iter->second);
FSET::const_iterator fset_iter;
for (fset_iter = fset->begin(); fset_iter != fset->end(); fset_iter++)
{
const PAIR &face = fset_iter->first;
const int elem = face.first;
// if this is not our element, do not do calculations
if (!feMesh_->is_owned_elt(elem)) continue;
const Array <XT_Function*> &fs = fset_iter->second;
_conn_ = feMesh_->element_connectivity_unique(elem);
// evaluate location at ips
feMesh_->face_shape_function(face, _fN_, _fdN_, _nN_, _fweights_);
feMesh_->element_coordinates(elem, xCoords);
MultAB(xCoords,_fN_,xAtIPs,0,1); //xAtIPs = xCoords*(N.transpose());
// interpolate prescribed flux at ips of this element
FSET::const_iterator face_iter = fset->find(face);
if (face_iter == fset->end())
{
stringstream ss;
ss << "face not found" << std::endl;
print_msg(communicator_,ss.str());
}
int nFieldDOF = (face_iter->second).size();
faceSource.reset(nIPsPerFace_,nFieldDOF);
for (int ip = 0; ip < nIPsPerFace_; ++ip) {
for (int idof = 0; idof<nFieldDOF; ++idof) {
XT_Function * f = fs(idof);
if (!f) continue;
faceSource(ip,idof) = f->f(column(xAtIPs,ip).ptr(),time);
}
}
// assemble
DENS_MAT & myNodalSources(nodalSources[_fieldName_].set_quantity());
_Nmat_ = _fN_.transMat(_fweights_*faceSource);
for (int i = 0; i < nNodesPerElement_; ++i)
{
int inode = _conn_(i);
for (int idof = 0; idof < nFieldDOF; ++idof)
{
myNodalSources(inode,idof) += _Nmat_(i,idof);
} // end assemble nFieldDOF
} // end assemble nNodesPerElement_
} // end fset loop
} // field loop
// assemble partial result matrices
for (src_iter=sourceFunctions.begin(); src_iter!=sourceFunctions.end(); src_iter++) {
_fieldName_ = src_iter->first;
if (!fieldMask((int)_fieldName_)) continue;
DENS_MAT & myNodalSources(nodalSources[_fieldName_].set_quantity());
allsum(communicator_,MPI_IN_PLACE, myNodalSources.ptr(), myNodalSources.size());
}
}
//-----------------------------------------------------------------
// prescribed volume flux using native quadrature
// integrate \int_V _N_I s(x,t) dV
//-----------------------------------------------------------------
void FE_Engine::add_sources(const Array<bool> &fieldMask,
const double time,
const VOLUME_SOURCE &sourceFunctions,
FIELDS &nodalSources) const
{
// sizing working arrays
DENS_MAT elemSource;
DENS_MAT xCoords(nSD_,nNodesPerElement_);
DENS_MAT xAtIPs(nSD_,nIPsPerElement_);
if (!(rank_zero(communicator_))) {
// Zero out unmasked nodal sources on all non-main processors.
// This is to avoid counting the previous nodal source values
// multiple times when we aggregate.
for (VOLUME_SOURCE::const_iterator src_iter = sourceFunctions.begin();
src_iter != sourceFunctions.end(); src_iter++) {
_fieldName_ = src_iter->first;
int index = (int) _fieldName_;
if (fieldMask(index)) {
DENS_MAT & myNodalSources(nodalSources[_fieldName_].set_quantity());
myNodalSources.reset(myNodalSources.nRows(), myNodalSources.nCols());
}
}
}
vector<int> myElems = feMesh_->owned_elts();
for (vector<int>::iterator elemsIter = myElems.begin();
elemsIter != myElems.end();
++elemsIter)
{
int ielem = *elemsIter;
_conn_ = feMesh_->element_connectivity_unique(ielem);
// evaluate location at ips
feMesh_->shape_function(ielem, _N_, _weights_);
feMesh_->element_coordinates(ielem, xCoords);
xAtIPs =xCoords*(_N_.transpose());
for (VOLUME_SOURCE::const_iterator src_iter = sourceFunctions.begin();
src_iter != sourceFunctions.end(); src_iter++) {
_fieldName_ = src_iter->first;
int index = (int) _fieldName_;
if ( fieldMask(index) ) {
const Array2D<XT_Function *> * thisSource
= (const Array2D<XT_Function *> *) &(src_iter->second);
int nFieldDOF = thisSource->nCols();
elemSource.reset(nIPsPerElement_,nFieldDOF);
// interpolate prescribed flux at ips of this element
for (int ip = 0; ip < nIPsPerElement_; ++ip) {
for (int idof = 0; idof < nFieldDOF; ++idof) {
XT_Function * f = (*thisSource)(ielem,idof);
if (f) {
elemSource(ip,idof) = f->f(column(xAtIPs,ip).ptr(),time);
}
}
}
// assemble
_Nmat_ = _N_.transMat(_weights_*elemSource);
DENS_MAT & myNodalSources(nodalSources[_fieldName_].set_quantity());
for (int i = 0; i < nNodesPerElement_; ++i) {
int inode = _conn_(i);
for (int idof = 0; idof < nFieldDOF; ++idof) {
myNodalSources(inode,idof) += _Nmat_(i,idof);
}
}
}
}
}
// Aggregate unmasked nodal sources on all processors.
for (VOLUME_SOURCE::const_iterator src_iter = sourceFunctions.begin();
src_iter != sourceFunctions.end(); src_iter++) {
_fieldName_ = src_iter->first;
int index = (int) _fieldName_;
if (fieldMask(index)) {
DENS_MAT & myNodalSources(nodalSources[_fieldName_].set_quantity());
allsum(communicator_,MPI_IN_PLACE, myNodalSources.ptr(), myNodalSources.size());
}
}
}
//-----------------------------------------------------------------
// boundary integral of a nodal field
//-----------------------------------------------------------------
void FE_Engine::field_surface_flux(
const DENS_MAT & field,
const set< PAIR > & faceSet,
DENS_MAT & values,
const bool contour,
const int axis) const
{
int dof = field.nCols();
double a[3] = {0,0,0};
a[axis] = 1;
// sizing working arrays
DENS_MAT n(nSD_,nIPsPerFace_);
DENS_MAT localElementFields(nNodesPerElement_,dof);
DENS_MAT integrals(dof,nSD_);
DENS_MAT fieldsAtIPs;
// SJL shouldn't this just be _fieldsAtIPs_
// the data member?
// element wise assembly
set< pair <int,int> >::iterator iter;
for (iter = faceSet.begin(); iter != faceSet.end(); iter++) {
int ielem = iter->first;
// if this is not our element, do not do calculations
//if (!feMesh_->is_owned_elt(ielem)) continue;
// evaluate shape functions at ips
feMesh_->face_shape_function(*iter, _N_, n, _fweights_);
// cross n with axis to get tangent
if (contour) {
double t[3];
for (int i = 0; i < nIPsPerFace_; i++) {
t[0] = a[1]*n(2,i) - a[2]*n(1,i);
t[1] = a[2]*n(0,i) - a[0]*n(2,i);
t[2] = a[0]*n(1,i) - a[1]*n(0,i);
n(0,i) = t[0];
n(1,i) = t[1];
n(2,i) = t[2];
}
}
// get connectivity
_conn_ = feMesh_->element_connectivity_unique(ielem);
// interpolate fields and gradients of fields ips of this element
for (int i = 0; i < nNodesPerElement_; i++) {
for (int j = 0; j < dof; j++) {
localElementFields(i,j) = field(_conn_(i),j);
}
}
// ips X dof = ips X nodes * nodes X dof
fieldsAtIPs = _N_*localElementFields;
// assemble : integral(k,j) = sum_ip n(j,ip) wg(ip,ip) field(ip,k)
_Nmat_ = n*_fweights_*fieldsAtIPs;
for (int j = 0; j < nSD_; j++) {
for (int k = 0; k < dof; ++k) {
integrals(k,j) += _Nmat_(j,k);
}
}
}
// assemble partial results
//LammpsInterface::instance()->allsum(MPI_IN_PLACE, integrals.ptr(), integrals.size());
// (S.n)_1 = S_1j n_j = S_11 n_1 + S_12 n_2 + S_13 n_3
// (S.n)_2 = S_2j n_j = S_21 n_1 + S_22 n_2 + S_23 n_3
// (S.n)_3 = S_3j n_j = S_31 n_1 + S_32 n_2 + S_33 n_3
if (dof == 9) { // tensor
values.reset(nSD_,1);
values(0,0) = integrals(0,0)+integrals(1,1)+integrals(2,2);
values(1,0) = integrals(3,0)+integrals(4,1)+integrals(5,2);
values(2,0) = integrals(6,0)+integrals(7,1)+integrals(8,2);
}
else if (dof == 6) { // sym tensor
values.reset(nSD_,1);
values(0,0) = integrals(0,0)+integrals(3,1)+integrals(4,2);
values(1,0) = integrals(3,0)+integrals(1,1)+integrals(5,2);
values(2,0) = integrals(4,1)+integrals(5,1)+integrals(2,2);
}
// (v.n) = v_j n_j
else if (dof == 3) { // vector
values.reset(1,1);
values(0,0) = integrals(0,0)+integrals(1,1)+integrals(2,2);
}
// s n = s n_j e_j
else if (dof == 1) { // scalar
values.reset(nSD_,1);
values(0,0) = integrals(0,0);
values(1,0) = integrals(0,1);
values(2,0) = integrals(0,2);
}
else {
string msg = "FE_Engine::field_surface_flux unsupported field width: ";
msg += to_string(dof);
throw ATC_Error(msg);
}
}
//-----------------------------------------------------------------
// evaluate shape functions at given points
//-----------------------------------------------------------------
void FE_Engine::evaluate_shape_functions(
const MATRIX & pt_coords,
SPAR_MAT & N) const
{
// Get shape function and derivatives at atomic locations
int nnodes = feMesh_->num_nodes_unique();
int npts = pt_coords.nRows();
// loop over point (atom) coordinates
DENS_VEC x(nSD_);
Array<int> node_index(nNodesPerElement_);
DENS_VEC shp(nNodesPerElement_);
N.reset(npts,nnodes);
for (int i = 0; i < npts; ++i) {
for (int k = 0; k < nSD_; ++k) { x(k) = pt_coords(i,k); }
feMesh_->shape_functions(x,shp,node_index);
for (int j = 0; j < nNodesPerElement_; ++j) {
int jnode = node_index(j);
N.add(i,jnode,shp(j));
}
}
N.compress();
}
//-----------------------------------------------------------------
// evaluate shape functions and their spatial derivatives at given points
//-----------------------------------------------------------------
void FE_Engine::evaluate_shape_functions(
const MATRIX & pt_coords,
SPAR_MAT & N,
SPAR_MAT_VEC & dN) const
{
// Get shape function and derivatives at atomic locations
int nnodes = feMesh_->num_nodes_unique();
int npts = pt_coords.nRows();
// loop over point (atom) coordinates
DENS_VEC x(nSD_);
Array<int> node_index(nNodesPerElement_);
DENS_VEC shp(nNodesPerElement_);
DENS_MAT dshp(nSD_,nNodesPerElement_);
for (int k = 0; k < nSD_; ++k) {
dN[k]->reset(npts,nnodes);
}
N.reset(npts,nnodes);
for (int i = 0; i < npts; ++i) {
for (int k = 0; k < nSD_; ++k) { x(k) = pt_coords(i,k); }
feMesh_->shape_functions(x,shp,dshp,node_index);
for (int j = 0; j < nNodesPerElement_; ++j) {
int jnode = node_index(j);
N.add(i,jnode,shp(j));
for (int k = 0; k < nSD_; ++k) {
dN[k]->add(i,jnode,dshp(k,j));
}
}
}
N.compress();
for (int k = 0; k < nSD_; ++k) {
dN[k]->compress();
}
}
//-----------------------------------------------------------------
// evaluate shape functions at given points
//-----------------------------------------------------------------
void FE_Engine::evaluate_shape_functions(
const MATRIX & pt_coords,
const INT_ARRAY & pointToEltMap,
SPAR_MAT & N) const
{
// Get shape function and derivatives at atomic locations
int nnodes = feMesh_->num_nodes_unique();
int npts = pt_coords.nRows();
// loop over point (atom) coordinates
DENS_VEC x(nSD_);
Array<int> node_index(nNodesPerElement_);
DENS_VEC shp(nNodesPerElement_);
N.reset(npts,nnodes);
for (int i = 0; i < npts; ++i) {
for (int k = 0; k < nSD_; ++k) { x(k) = pt_coords(i,k); }
int eltId = pointToEltMap(i,0);
feMesh_->shape_functions(x,eltId,shp,node_index);
for (int j = 0; j < nNodesPerElement_; ++j) {
int jnode = node_index(j);
N.add(i,jnode,shp(j));
}
}
N.compress();
}
//-----------------------------------------------------------------
// evaluate shape functions and their spatial derivatives at given points
//-----------------------------------------------------------------
void FE_Engine::evaluate_shape_functions(
const MATRIX & pt_coords,
const INT_ARRAY & pointToEltMap,
SPAR_MAT & N,
SPAR_MAT_VEC & dN) const
{
// Get shape function and derivatives at atomic locations
int nnodes = feMesh_->num_nodes_unique();
int npts = pt_coords.nRows();
// loop over point (atom) coordinates
DENS_VEC x(nSD_);
Array<int> node_index(nNodesPerElement_);
DENS_VEC shp(nNodesPerElement_);
DENS_MAT dshp(nSD_,nNodesPerElement_);
for (int k = 0; k < nSD_; ++k) {
dN[k]->reset(npts,nnodes);
}
N.reset(npts,nnodes);
for (int i = 0; i < npts; ++i) {
for (int k = 0; k < nSD_; ++k) { x(k) = pt_coords(i,k); }
int eltId = pointToEltMap(i,0);
feMesh_->shape_functions(x,eltId,shp,dshp,node_index);
for (int j = 0; j < nNodesPerElement_; ++j) {
int jnode = node_index(j);
N.add(i,jnode,shp(j));
for (int k = 0; k < nSD_; ++k) {
dN[k]->add(i,jnode,dshp(k,j));
}
}
}
N.compress();
for (int k = 0; k < nSD_; ++k) {
dN[k]->compress();
}
}
//-----------------------------------------------------------------
// evaluate shape functions spatial derivatives at given points
//-----------------------------------------------------------------
void FE_Engine::evaluate_shape_function_derivatives(
const MATRIX & pt_coords,
const INT_ARRAY & pointToEltMap,
SPAR_MAT_VEC & dNdx ) const
{
int nnodes = feMesh_->num_nodes_unique();
int npts = pt_coords.nRows();
for (int k = 0; k < nSD_; ++k) {
dNdx[k]->reset(npts,nnodes);
}
// loop over point (atom) coordinates
DENS_VEC x(nSD_);
Array<int> node_index(nNodesPerElement_);
DENS_MAT dshp(nSD_,nNodesPerElement_);
for (int i = 0; i < npts; ++i) {
for (int k = 0; k < nSD_; ++k) { x(k) = pt_coords(i,k); }
int eltId = pointToEltMap(i,0);
feMesh_->shape_function_derivatives(x,eltId,dshp,node_index);
for (int j = 0; j < nNodesPerElement_; ++j) {
int jnode = node_index(j);
for (int k = 0; k < nSD_; ++k) {
dNdx[k]->add(i,jnode,dshp(k,j));
}
}
}
for (int k = 0; k < nSD_; ++k) {
dNdx[k]->compress();
}
}
//-------------------------------------------------------------------
// set kernels
//-------------------------------------------------------------------
void FE_Engine::set_kernel(KernelFunction* ptr){kernelFunction_=ptr;}
//-------------------------------------------------------------------
// kernel_matrix_bandwidth
//-------------------------------------------------------------------
int FE_Engine::kernel_matrix_bandwidth(
const MATRIX & ptCoords) const
{
if (! kernelFunction_) throw ATC_Error("no kernel function specified");
int npts = ptCoords.nRows();
int bw = 0;
if (npts>0) {
DENS_VEC xI(nSD_),x(nSD_),dx(nSD_);
for (int i = 0; i < nNodes_; i++) {
xI = feMesh_->global_coordinates(i);
for (int j = 0; j < npts; j++) {
for (int k = 0; k < nSD_; ++k) { x(k) = ptCoords(j,k); }
dx = x - xI;
if (kernelFunction_->in_support(dx)) bw = max(bw,abs(j-map_global_to_unique(i)));
}
}
}
return bw;
}
//-------------------------------------------------------------------
// evaluate kernels
//-------------------------------------------------------------------
void FE_Engine::evaluate_kernel_functions(
const MATRIX & ptCoords,
SPAR_MAT & N ) const
{
print_msg_once(communicator_,"kernel matrix bandwidth " +to_string(kernel_matrix_bandwidth(ptCoords)));
if (! kernelFunction_) throw ATC_Error("no kernel function specified");
int npts = ptCoords.nRows();
if (npts>0) {
N.reset(npts,nNodesUnique_); // transpose of N_Ia
DENS_VEC xI(nSD_),x(nSD_),dx(nSD_);
for (int i = 0; i < nNodes_; i++) {
xI = feMesh_->global_coordinates(i);
for (int j = 0; j < npts; j++) {
for (int k = 0; k < nSD_; ++k) { x(k) = ptCoords(j,k); }
dx = x - xI;
double val = kernelFunction_->value(dx);
val *= kernelFunction_->dimensionality_factor();
if (val > 0) N.add(j,map_global_to_unique(i),val);
}
}
N.compress();
}
}
//-----------------------------------------------------------------
// create a non-uniform rectangular mesh on a rectangular region
//-----------------------------------------------------------------
void FE_Engine::create_mesh(Array<double> & dx,
Array<double> & dy,
Array<double> & dz,
const char * regionName,
Array<bool> periodicity)
{
double xmin, xmax, ymin, ymax, zmin, zmax;
double xscale, yscale, zscale;
// check to see if region exists and is it a box, and if so get the bounds
bool isBox;
isBox = ATC::LammpsInterface::instance()->region_bounds(
regionName,
xmin, xmax,
ymin, ymax,
zmin, zmax,
xscale,
yscale,
zscale);
if (!isBox) throw ATC_Error("Region for FE mesh is not a box");
if (dx(0) == 0.0) dx = (xmax-xmin)/dx.size();
if (dy(0) == 0.0) dy = (ymax-ymin)/dy.size();
if (dz(0) == 0.0) dz = (zmax-zmin)/dz.size();
feMesh_ = new FE_Rectangular3DMesh(dx,dy,dz,
xmin, xmax,
ymin, ymax,
zmin, zmax,
periodicity,
xscale, yscale, zscale);
stringstream ss;
ss << "created structured mesh with " << feMesh_->num_nodes() << " nodes, "
<< feMesh_->num_nodes_unique() << " unique nodes, and "
<< feMesh_->num_elements() << " elements";
print_msg_once(communicator_,ss.str());
#ifdef ATC_VERBOSE
string msg = "element sizes in x:\n";
double sum = 0;
for (int i = 0; i < dx.size(); ++i) { sum += dx(i); }
double xn= 1.0/sum;
double xs= xn*(xmax-xmin);
double xl= xs/(ATC::LammpsInterface::instance()->xlattice());
double sumn = 0, sums = 0, suml = 0;
for (int i = 0; i < dx.size(); ++i) {
double dxn = dx(i)*xn; sumn += dxn;
double dxs = dx(i)*xs; sums += dxs;
double dxl = dx(i)*xl; suml += dxl;
msg += " "+to_string(i+1)
+" "+to_string(dx(i))
+" "+to_string(dxn)
+" "+to_string(dxs)
+" "+to_string(dxl)+"\n";
}
msg += "sum "+to_string(sum)+" "
+to_string(sumn)+" "
+to_string(sums)+" "
+to_string(suml)+"\n";
print_msg_once(communicator_,msg);
#endif
}
//-----------------------------------------------------------------
// create a uniform rectangular mesh on a rectangular region
//-----------------------------------------------------------------
void FE_Engine::create_mesh(int nx, int ny, int nz, const char * regionName,
Array<bool> periodicity)
{
double xmin, xmax, ymin, ymax, zmin, zmax;
double xscale, yscale, zscale;
// check to see if region exists and is it a box, and if so get the bounds
bool isBox;
isBox = ATC::LammpsInterface::instance()->region_bounds(
regionName,
xmin, xmax,
ymin, ymax,
zmin, zmax,
xscale, yscale, zscale);
if (!isBox) throw ATC_Error("Region for FE mesh is not a box");
feMesh_ = new FE_Uniform3DMesh(nx,ny,nz,
xmin, xmax,
ymin, ymax,
zmin, zmax,
periodicity,
xscale, yscale, zscale);
stringstream ss;
ss << "created uniform mesh with " << feMesh_->num_nodes() << " nodes, "
<< feMesh_->num_nodes_unique() << " unique nodes, and "
<< feMesh_->num_elements() << " elements";
print_msg_once(communicator_,ss.str());
}
//-----------------------------------------------------------------
// read a mesh from an input file using MeshReader object
//-----------------------------------------------------------------
void FE_Engine::read_mesh(string meshFile, Array<bool> & periodicity)
{
if (feMesh_) throw ATC_Error("FE_Engine already has a mesh");
feMesh_ = MeshReader(meshFile,periodicity).create_mesh();
feMesh_->initialize();
}
};
Event Timeline
Log In to Comment