Homec4science
Diffusion lammps (lammm-master) Browse

lammps/lib/atc/FE_Engine.cpp
dfb38392c991lammm-master

Last Change
Raw File

FE_Engine.cpp
View Options

#include "FE_Engine.h"
#include "ATC_Transfer.h"
#include "FE_Element.h"
#include "XT_Function.h"
#include "LammpsInterface.h"
#include "PhysicsModel.h"
#include <stdio.h>
#include <map>
using namespace std;
namespace ATC{
//-----------------------------------------------------------------
FE_Engine::FE_Engine(ATC_Transfer * atcTransfer)
: atcTransfer_(atcTransfer),
feMesh_(NULL),
initialized_(false),
amendedMeshData_(false),
outputManager_()
{
}
FE_Engine::~FE_Engine()
{
if (feMesh_) delete feMesh_;
if (connectivity_ && amendedMeshData_) delete connectivity_;
if (nodeMap_ && amendedMeshData_) delete nodeMap_;
if (coordinates_ && amendedMeshData_) delete coordinates_;
}
//-----------------------------------------------------------------
void FE_Engine::initialize()
{
if (! feMesh_)
throw ATC_Error(0,"FE_Engine has no mesh");
if (! initialized_ ) {
// mesh data
nodeMap_ = &(feMesh_->global_to_unique_map());
connectivity_ = &(feMesh_->connectivity());
coordinates_ = &(feMesh_->nodal_coordinates());
// set up work spaces
nNodesPerElement_ = feMesh_->get_nNodesPerElement();
nIPsPerElement_ = feMesh_->get_nIPsPerElement();
nIPsPerFace_ = feMesh_->get_nIPsPerFace();
nSD_ = feMesh_->get_nSpatialDimensions();
nElems_ = feMesh_->get_nElements();
// remove specified elements
if (nullElements_.size() > 0) delete_elements(nullElements_);
initialized_ = true;
}
}
//-----------------------------------------------------------------
bool FE_Engine::modify(int narg, char **arg)
{
bool match = false;
/*! \page man_fem_mesh fix_modify AtC fem create mesh
\section syntax
fix_modify AtC fem create mesh <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
\section examples
<TT> fix_modify AtC fem 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
*/
if (strcmp(arg[0],"fem")==0) {
// create mesh
if (strcmp(arg[1],"create")==0) {
if (strcmp(arg[2],"mesh")==0) {
int nx = atoi(arg[3]);
int ny = atoi(arg[4]);
int nz = atoi(arg[5]);
int xperiodic = 0;
if (strcmp(arg[7],"p")==0) xperiodic = 1;
int yperiodic = 0;
if (strcmp(arg[8],"p")==0) yperiodic = 1;
int zperiodic = 0;
if (strcmp(arg[9],"p")==0) zperiodic = 1;
if (feMesh_) throw ATC_Error(0,"FE Engine already has a mesh");
create_mesh(nx, ny, nz, arg[6], xperiodic, yperiodic, zperiodic );
// construct prescribed data manager NOTE move to ATC_Transfer?
atcTransfer_->initialize_mesh_data();
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[0],"mesh")==0)
&& (strcmp(arg[1],"delete_elements")==0) ) {
string esetName = arg[2];
set<int> elemSet = feMesh_->get_elementset(esetName);
nullElements_.insert(elemSet.begin(), elemSet.end());
match = true;
}
// FE_Mesh
else {
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, OutputType otype)
{
outputManager_.initialize(outputPrefix, otype);
if (! feMesh_) throw ATC_Error(0,"output needs mesh");
if (! initialized_) initialize();
if (rank == 0) outputManager_.write_geometry(*coordinates_, connectivity_);
}
//-----------------------------------------------------------------
// write geometry
//-----------------------------------------------------------------
void FE_Engine::write_geometry(void)
{
outputManager_.write_geometry(*coordinates_, connectivity_);
}
// -------------------------------------------------------------
// write data
// -------------------------------------------------------------
void FE_Engine::write_data(double time, FIELDS &soln, OUTPUT_LIST *data)
{
outputManager_.write_data(time, &soln, data, nodeMap_->get_data());
}
// -------------------------------------------------------------
// write data
// -------------------------------------------------------------
void FE_Engine::write_data(double time, OUTPUT_LIST *data)
{
outputManager_.write_data(time, data, nodeMap_->get_data());
}
// -------------------------------------------------------------
// amend mesh for deleted elements
// -------------------------------------------------------------
void FE_Engine::delete_elements(const set<int> & elementList)
{
int nsd = feMesh_->get_nSpatialDimensions();
set<int> elementsNew;
feMesh_->elementset_complement(elementList,elementsNew);
int nElementsNew = elementsNew.size();
set<int> newToOld;
map<int,int> oldToNewMap;
feMesh_->elementset_to_nodeset(elementsNew,newToOld);
int nNodesNew = newToOld.size();
set<int>::const_iterator itr;
// coordinates & node map (from nodes to data)
const DENS_MAT &coordinates = feMesh_->nodal_coordinates();
const Array<int> & node_map = feMesh_->global_to_unique_map();
if (coordinates_ && amendedMeshData_) delete coordinates_;
coordinates_ = new DENS_MAT(nsd,nNodesNew);
DENS_MAT & coor = * (const_cast<DENS_MAT *> ( coordinates_));
if (nodeMap_ && amendedMeshData_) delete nodeMap_;
nodeMap_ = new Array<int> (nNodesNew);
Array<int> & nmap = * (const_cast<Array<int> *> ( nodeMap_));
int k = 0, i = 0;
for (itr = newToOld.begin(); itr != newToOld.end(); itr++) {
int node = *itr;
oldToNewMap[node]=i++;
nmap(k) = node_map(node);
for(int j = 0; j < nsd; j++) {
coor(j,k) = coordinates(j,node);
}
k++;
}
// connectivity
const Array2D<int> & connectivity = feMesh_->connectivity();
int nNodesPerElement = connectivity.nRows();
int nElements = connectivity.nCols();
if (connectivity_ && amendedMeshData_) delete connectivity_;
connectivity_ = new Array2D<int>(nNodesPerElement,nElementsNew);
Array2D<int> & conn = * (const_cast<Array2D<int> *> ( connectivity_));
k = 0;
for (itr = elementsNew.begin(); itr != elementsNew.end(); itr++) {
int ielem = *itr;
for(int j = 0; j < nNodesPerElement; j++) {
int old_node = connectivity(j,ielem);
map<int,int>::iterator map_itr = oldToNewMap.find(old_node);
if (map_itr == oldToNewMap.end()) {
cout << "map failure " << old_node << "\n";
}
int node = map_itr->second;
conn(j,k) = node;
}
k++;
}
amendedMeshData_ = true;
}
// -------------------------------------------------------------
// compute dimensionless stiffness matrix using native quadrature
// -------------------------------------------------------------
void FE_Engine::compute_stiffness_matrix(SPAR_MAT &matrix) const
{
//int nfields = field_mask.get_length();
int nNodes = feMesh_->get_nNodesUnique();
int nNodesPerElement = feMesh_->get_nNodesPerElement();
int nIPsPerElement = feMesh_->get_nIPsPerElement();
int nsd = feMesh_->get_nSpatialDimensions();
int nelems = feMesh_->get_nElements();
DENS_MAT N(nIPsPerElement,nNodesPerElement);
vector<DENS_MAT> dN(nsd);
dN.assign(nsd, DENS_MAT(nIPsPerElement,nNodesPerElement));
DiagonalMatrix<double> weights(nIPsPerElement,nIPsPerElement);
Array<int> conn(nNodesPerElement);
// assemble consistent mass (nnodes X nnodes)
matrix.reset(nNodes,nNodes);
DENS_MAT elementMassMatrix(nNodesPerElement,nNodesPerElement);
for (int ielem = 0; ielem < nelems; ++ielem)
{
// evaluate shape functions
feMesh_->shape_function(ielem, N, dN, weights);
// 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
feMesh_->element_connectivity_unique(ielem, conn);
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));
}
}
}
}
// -------------------------------------------------------------
// compute dimensionless consistent mass using native quadrature
// -------------------------------------------------------------
void FE_Engine::compute_mass_matrix(SPAR_MAT &mass_matrix) const
{
int nNodes = feMesh_->get_nNodesUnique();
int nNodesPerElement = feMesh_->get_nNodesPerElement();
int nIPsPerElement = feMesh_->get_nIPsPerElement();
int nsd = feMesh_->get_nSpatialDimensions();
int nelems = feMesh_->get_nElements();
DENS_MAT N(nIPsPerElement,nNodesPerElement);
DiagonalMatrix<double> weights(nIPsPerElement,nIPsPerElement);
Array<int> conn(nNodesPerElement);
// assemble consistent mass (nnodes X nnodes)
mass_matrix.reset(nNodes,nNodes);
DENS_MAT elementMassMatrix(nNodesPerElement,nNodesPerElement);
for (int ielem = 0; ielem < nelems; ++ielem)
{
// evaluate shape functions
feMesh_->shape_function(ielem, N, weights);
// perform quadrature
elementMassMatrix = N.transMat(weights*N);
// get connectivity
feMesh_->element_connectivity_unique(ielem, conn);
for (int i = 0; i < nNodesPerElement; ++i)
{
int inode = conn(i);
for (int j = 0; j < nNodesPerElement; ++j)
{
int jnode = conn(j);
mass_matrix.add(inode, jnode, elementMassMatrix(i,j));
}
}
}
}
// -------------------------------------------------------------
// compute consistent mass using given quadrature
// -------------------------------------------------------------
void FE_Engine::compute_mass_matrix(const DIAG_MAT &weights,
const SPAR_MAT &N,
SPAR_MAT &mass_matrix) 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;
LammpsInterface::instance()->allsum(
tmp_mass_matrix_local.get_ptr(),
tmp_mass_matrix.get_ptr(), count);
// create sparse from dense
mass_matrix.reset(tmp_mass_matrix);
}
// -------------------------------------------------------------
// compute dimensionless lumped mass using native quadrature
// -------------------------------------------------------------
void FE_Engine::compute_lumped_mass_matrix(DiagonalMatrix<double> & M) const
{
// initialize
int nNodes = feMesh_->get_nNodesUnique();
M.reset(nNodes,nNodes);
int nNodesPerElement = feMesh_->get_nNodesPerElement();
int nIPsPerElement = feMesh_->get_nIPsPerElement();
int nsd = feMesh_->get_nSpatialDimensions();
int nelems = feMesh_->get_nElements();
DENS_MAT N(nNodesPerElement,nIPsPerElement);
DIAG_MAT weights(nIPsPerElement,nIPsPerElement);
Array<int> conn(nNodesPerElement);
// assemble lumped diagonal mass (nnodes X nnodes)
int inode = -1;
DENS_VEC weightVec(nIPsPerElement);
DENS_VEC Nvec(nNodesPerElement);
for (int ielem = 0; ielem < nelems; ++ielem)
{
// evaluate shape functions
feMesh_->shape_function(ielem, N, weights);
int nips = weights.nRows(); // different than nIPsPerElement?
weightVec = 1.;
weightVec = weights*weightVec;
// get connectivity
feMesh_->element_connectivity_unique(ielem, conn);
Nvec = N*weightVec;
for (int i = 0; i < nNodesPerElement; ++i)
{
inode = conn(i);
M(inode,inode) += Nvec(i);
}
}
}
// -------------------------------------------------------------
// 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,
map<FieldName, DIAG_MAT>& M, // mass matrix
const Array<bool> *element_mask) const
{
int nfields = field_mask.get_length();
// initialize
int nNodes = feMesh_->get_nNodesUnique();
for (int j = 0; j < nfields; ++j) {
M[field_mask(j)].reset(nNodes,nNodes);
}
int nNodesPerElement = feMesh_->get_nNodesPerElement();
int nIPsPerElement = feMesh_->get_nIPsPerElement();
int nsd = feMesh_->get_nSpatialDimensions();
int nelems = feMesh_->get_nElements();
DENS_MAT N(nNodesPerElement,nIPsPerElement);
DIAG_MAT weights(nIPsPerElement,nIPsPerElement);
Array<int> conn(nNodesPerElement);
FIELDS fieldsAtIPs, localElementFields;
DENS_MAT Nmat;
// assemble lumped diagonal mass (nnodes X nnodes)
for (int ielem = 0; ielem < nelems; ++ielem)
{
// if element is masked, skip it
if (element_mask && !(*element_mask)(ielem)) continue;
// material id
int imat = elementMaterials(ielem);
// evaluate shape functions
feMesh_->shape_function(ielem, N, weights);
int nips = weights.nRows();
// get connectivity
feMesh_->element_connectivity_unique(ielem, conn);
// interpolate fields at IPs
FIELDS::const_iterator field;
for (field = fields.begin(); field != fields.end(); field++)
{
// field values at all nodes
const DENS_MAT &vI = field->second;
// field values at element nodes
DENS_MAT &vIe = localElementFields[field->first];
// field values at integration points -> to be computed
DENS_MAT &vP = fieldsAtIPs[field->first];
int numFieldDOF = vI.nCols();
vIe.reset(nNodesPerElement, numFieldDOF);
// gather local field
for (int i = 0; i < nNodesPerElement; i++) {
for (int j = 0; j < numFieldDOF; j++) {
vIe(i,j) = vI(conn(i),j);
}
}
// interpolate field at integration points
vP = N*vIe;
}
// compute energy/momentum/mass densities
FIELDS capacities;
physicsModel->M_integrand(field_mask, fieldsAtIPs, capacities, imat);
// integrate & assemble
for (int j = 0; j < nfields; ++j) {
FieldName thisFieldName = field_mask(j);
Nmat = N.transMat(weights*capacities[thisFieldName]);
for (int i = 0; i < nNodesPerElement; ++i) {
int inode = conn(i);
M[thisFieldName](inode,inode) += Nmat(i,0);// assume all dof same
}
}
}
}
// -------------------------------------------------------------
// 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, // NOTE use these as a mask
const SPAR_MAT &N,
map<FieldName, DIAG_MAT> &M) const // mass matrices
{
int nips = weights.nCols();
int nsd = feMesh_->get_nSpatialDimensions();
int nfields = field_mask.get_length();
int nNodes = feMesh_->get_nNodesUnique();
// initialize
FieldName thisFieldName;
map<FieldName, DIAG_MAT> M_local;
for (int j = 0; j < nfields; ++j) {
thisFieldName = field_mask(j);
M_local[thisFieldName].reset(nNodes,nNodes);
M[thisFieldName].reset(nNodes,nNodes);
}
if (nips>0) {
// compute fields at given ips
FIELDS fieldsAtIPs;
for (int j = 0; j < nfields; ++j) {
thisFieldName = field_mask(j);
const FIELD & field = (fields.find(thisFieldName))->second;
int numFieldDOF = field.nCols();
fieldsAtIPs[thisFieldName].reset(nips,numFieldDOF);
fieldsAtIPs[thisFieldName] = 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++;
}
bool singleMaterial = ( atomMatls == 1 );
// setup data structures
FIELDS capacities;
// evaluate physics model
if (singleMaterial) {
physicsModel->M_integrand(field_mask, fieldsAtIPs, capacities);
}
else {
throw ATC_Error(0,"unimplemented function in FE_Engine::compute_mass_matrix");
}
// integrate & assemble
for (int j = 0; j < nfields; ++j) {
FieldName thisFieldName = field_mask(j);
M_local[thisFieldName].reset( // assume all columns same
column(N.transMat(weights*capacities[thisFieldName]),0) );
}
}
// Share information between processors
for (int j = 0; j < nfields; ++j) {
FieldName thisFieldName = field_mask(j);
int count = M_local[thisFieldName].size();
LammpsInterface::instance()->allsum(
M_local[thisFieldName].get_ptr(),
M[thisFieldName].get_ptr(), count);
}
}
//-----------------------------------------------------------------
// compute assembled average gradient evaluated at the nodes
//-----------------------------------------------------------------
void FE_Engine::compute_gradient_matrix(
GRAD_SHPFCN & grad_matrix) const
{
int nNodesUnique = feMesh_->get_nNodesUnique();
int nNodesPerElement = feMesh_->get_nNodesPerElement();
int nsd = feMesh_->get_nSpatialDimensions();
int nelems = feMesh_->get_nElements();
int nIPsPerElement = feMesh_->get_nIPsPerElement();
DENS_MAT N(nIPsPerElement,nNodesPerElement);
vector< DENS_MAT > dN(nsd);
vector< DENS_MAT > tmp_grad_matrix(nsd);
for (int i = 0; i < nsd; i++) {
dN[i].reset(nNodesPerElement,nNodesPerElement);
tmp_grad_matrix[i].reset(nNodesUnique,nNodesUnique);
}
DiagonalMatrix<double> weights(nIPsPerElement,nIPsPerElement);
Array<int> conn(nNodesPerElement);
// element wise assembly
Array<int> count(nNodesUnique); count = 0;
feMesh_->set_quadrature(FE_Element::NODAL_QUADRATURE);
for (int ielem = 0; ielem < nelems; ++ielem) {
// evaluate shape functions
feMesh_->shape_function(ielem, N, dN, weights);
int nips = weights.nRows();
// get connectivity
feMesh_->element_connectivity_unique(ielem, conn);
for (int j = 0; j < nIPsPerElement; ++j) {
int jnode = conn(j); // NOTE assumes ip order is node order
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);
}
}
}
}
feMesh_->set_quadrature(FE_Element::GAUSSIAN_QUADRATURE); // reset to default
for (int inode = 0; inode < nNodesUnique; ++inode) {
//cout << inode << " count " << count(inode) << "\n";
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
grad_matrix.resize(nsd);
for (int k = 0; k < nsd; ++k) {
grad_matrix[k].reset(tmp_grad_matrix[k]);
}
}
// -------------------------------------------------------------
// 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,
FIELDS &energy,
const Array<bool> *element_mask) const
{
//NOTE move to initialization and make workspace
const int nNodesPerElement = feMesh_->get_nNodesPerElement();
const int nIPsPerElement = feMesh_->get_nIPsPerElement();
const int nsd = feMesh_->get_nSpatialDimensions();
const int nelems = feMesh_->get_nElements();
DENS_MAT Nmat;
//NOTE move to initialization and make workspace
DENS_MAT N(nIPsPerElement,nNodesPerElement);
vector<DENS_MAT> dN(nsd);
dN.assign(nsd, DENS_MAT(nIPsPerElement,nNodesPerElement));
FIELDS fieldsAtIPs, localElementFields;
GRAD_FIELDS grad_fieldsAtIPs;
FIELDS::const_iterator field;
FieldName thisFieldName;
int numFieldDOF;
for (int n = 0; n < mask.get_length(); n++) {
FieldName thisFieldName = mask(n);
const FIELD & field = (fields.find(thisFieldName))->second;
energy[thisFieldName].reset(field.nRows(), field.nCols());
}
DIAG_MAT weights(nIPsPerElement,nIPsPerElement);
Array<int> conn(nNodesPerElement);
// element wise assembly
int inode = -1;
for (int ielem=0; ielem < nelems; ielem++)
{
// if element is masked, skip it
if (element_mask && !(*element_mask)(ielem)) continue;
// material id
int imat = elementMaterials(ielem);
// evaluate shape functions
feMesh_->shape_function(ielem, N, dN, weights);
int nips = weights.nRows();
// get connectivity
feMesh_->element_connectivity_unique(ielem, conn);
// compute fields and gradients of fields ips of this element
for (int n = 0; n < mask.get_length(); n++) {
FieldName thisFieldName = mask(n);
FIELDS::const_iterator field = (fields.find(thisFieldName));
// field values at all nodes
const DENS_MAT &vI = field->second;
// field values at element nodes
DENS_MAT &vIe = localElementFields[field->first];
// field values at integration points -> to be computed
DENS_MAT &vP = fieldsAtIPs[field->first];
// gradients of field at integration points -> to be computed
vector<DENS_MAT> &dvP = grad_fieldsAtIPs[field->first];
numFieldDOF = vI.nCols();
vIe.reset(nNodesPerElement, numFieldDOF);
// gather local field
for (int i = 0; i < nNodesPerElement; i++)
for (int j = 0; j < numFieldDOF; j++)
vIe(i,j) = vI(conn(i),j);
// interpolate field at integration points
vP = N*vIe;
dvP.assign(nsd, DENS_MAT(nIPsPerElement, numFieldDOF));
for (int j = 0; j < nsd; ++j) dvP[j] = dN[j]*vIe;
}
FIELDS energies;
physicsModel->E_integrand(mask, fieldsAtIPs, grad_fieldsAtIPs, energies, imat);
// assemble
for (int n = 0; n < mask.get_length(); n++) {
FieldName thisFieldName = mask(n);
FIELDS::const_iterator field = (fields.find(thisFieldName));
numFieldDOF = (field->second).nCols();
Nmat = N.transMat(weights*energies[thisFieldName]);
for (int i = 0; i < nNodesPerElement; ++i) {
inode = conn(i);
for (int k = 0; k < numFieldDOF; ++k) {
energy[thisFieldName](inode,k) += Nmat(i,k);
}
}
}
} // element loop
} // function call
// -------------------------------------------------------------
// compute rhs using native quadrature
// -------------------------------------------------------------
void FE_Engine::compute_rhs_vector(const Array2D<bool> &rhs_mask,
const FIELDS &fields,
const PhysicsModel * physicsModel,
const Array<int> & elementMaterials,
FIELDS &rhs,
const Array<bool> *element_mask) const
{
//NOTE move to initialization and make workspace
const int nNodesPerElement = feMesh_->get_nNodesPerElement();
const int nIPsPerElement = feMesh_->get_nIPsPerElement();
const int nsd = feMesh_->get_nSpatialDimensions();
const int nelems = feMesh_->get_nElements();
DENS_MAT Nmat;
//NOTE move to initialization and make workspace
DENS_MAT N(nIPsPerElement,nNodesPerElement);
vector<DENS_MAT> dN(nsd);
dN.assign(nsd, DENS_MAT(nIPsPerElement,nNodesPerElement));
FIELDS fieldsAtIPs, localElementFields;
GRAD_FIELDS grad_fieldsAtIPs;
FIELDS::const_iterator field;
FieldName thisFieldName;
int numFieldDOF;
for (field = fields.begin(); field != fields.end(); field++) {
thisFieldName = field->first;
if (rhs_mask(thisFieldName,FLUX) || rhs_mask(thisFieldName,SOURCE)) {
int nrows = (field->second).nRows();
int ncols = (field->second).nCols();
rhs[thisFieldName].reset(nrows,ncols);
}
}
DIAG_MAT weights(nIPsPerElement,nIPsPerElement);
Array<int> conn(nNodesPerElement);
// element wise assembly
int inode = -1;
for (int ielem=0; ielem < nelems; ielem++)
{
// if element is masked, skip it
if (element_mask && !(*element_mask)(ielem)) continue;
// material id
int imat = elementMaterials(ielem);
// evaluate shape functions
feMesh_->shape_function(ielem, N, dN, weights);
int nips = weights.nRows();
// get connectivity
feMesh_->element_connectivity_unique(ielem, conn);
// compute fields and gradients of fields ips of this element
for (field = fields.begin(); field != fields.end(); field++)
{
// field values at all nodes
const DENS_MAT &vI = field->second;
// field values at element nodes
DENS_MAT &vIe = localElementFields[field->first];
// field values at integration points -> to be computed
DENS_MAT &vP = fieldsAtIPs[field->first];
// gradients of field at integration points -> to be computed
vector<DENS_MAT> &dvP = grad_fieldsAtIPs[field->first];
numFieldDOF = vI.nCols();
vIe.reset(nNodesPerElement, numFieldDOF);
// gather local field
for (int i = 0; i < nNodesPerElement; i++)
for (int j = 0; j < numFieldDOF; j++)
vIe(i,j) = vI(conn(i),j);
// interpolate field at integration points
vP = N*vIe;
dvP.assign(nsd, DENS_MAT(nIPsPerElement, numFieldDOF));
for (int j = 0; j < nsd; ++j) dvP[j] = dN[j]*vIe;
}
// NOTE move scaling of N by weights up here
// evaluate physics model
// N_fluxes is a set of fields
if(physicsModel->has_N_integrand()) {
FIELDS N_fluxes;
physicsModel->N_integrand(rhs_mask, fieldsAtIPs, grad_fieldsAtIPs, N_fluxes, imat);
// assemble
for (field = fields.begin(); field != fields.end(); field++)
{
thisFieldName = field->first;
// NOTE why is this line here?
// SOURCE should refer only to adding extrinsic coupling sources
if (!rhs_mask(thisFieldName,SOURCE)) continue;
numFieldDOF = (field->second).nCols();
Nmat = N.transMat(weights*N_fluxes[thisFieldName]);
for (int i = 0; i < nNodesPerElement; ++i) {
inode = conn(i);
for (int k = 0; k < numFieldDOF; ++k) {
rhs[thisFieldName](inode,k) += Nmat(i,k);
}
}
}
}
// evaluate Physics model
// B_fluxes is a set of field gradients
if (physicsModel->has_B_integrand()) {
GRAD_FIELDS B_fluxes;
physicsModel->B_integrand(rhs_mask, fieldsAtIPs, grad_fieldsAtIPs, B_fluxes, imat);
// assemble
for (field = fields.begin(); field != fields.end(); field++)
{
thisFieldName = field->first;
numFieldDOF = (field->second).nCols();
if (!rhs_mask(thisFieldName,FLUX)) continue;
for (int j = 0; j < nsd; j++)
{
Nmat = dN[j].transMat(weights*B_fluxes[thisFieldName][j]);
for (int i = 0; i < nNodesPerElement; ++i)
{
inode = conn(i);
for (int k = 0; k < numFieldDOF; ++k)
{
rhs[thisFieldName](inode,k) += Nmat(i,k);
} // scatter k
} // scatter i
} // loop over nsd
} //loop on fields
} // if B_integrand
} // element loop
} // function call
// -------------------------------------------------------------
// compute rhs using given quadrature
// -------------------------------------------------------------
void FE_Engine::compute_rhs_vector(const Array2D<bool> &rhs_mask,
const FIELDS &fields,
const PhysicsModel *physicsModel,
const Array<set<int> >&pointMaterialGroups,
const DIAG_MAT &weights,
const SPAR_MAT &N,
const GRAD_SHPFCN &dN,
FIELDS &rhs) const
{
int nips = weights.nCols();
int nNodesUnique = feMesh_->get_nNodesUnique();
int nsd = feMesh_->get_nSpatialDimensions();
int nelems = feMesh_->get_nElements();
FieldName thisFieldName;
FIELDS::const_iterator field;
FIELDS rhs_local;
for (field = rhs.begin(); field != rhs.end(); field++) {
thisFieldName = field->first;
if (rhs_mask(thisFieldName,FLUX) || rhs_mask(thisFieldName,SOURCE)) {
int nrows = (field->second).nRows();
int ncols = (field->second).nCols();
rhs [thisFieldName].reset(nrows,ncols);
rhs_local[thisFieldName].reset(nrows,ncols);
}
}
if (nips>0) {
// compute fields and gradients of fields at given ips
GRAD_FIELDS grad_fieldsAtIPs;
FIELDS fieldsAtIPs;
int numFieldDOF;
for (field = fields.begin(); field != fields.end(); field++) {
thisFieldName = field->first;
const FIELD & field = (fields.find(thisFieldName))->second;
numFieldDOF = field.nCols();
grad_fieldsAtIPs[thisFieldName].assign(nsd,DENS_MAT(nips,numFieldDOF));
fieldsAtIPs[thisFieldName].reset(nips,numFieldDOF);
fieldsAtIPs[thisFieldName] = N*field;
for (int j = 0; j < nsd; ++j) {
grad_fieldsAtIPs[thisFieldName][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 );
// setup data structures
FIELDS N_fluxes;
GRAD_FIELDS B_fluxes;
if(physicsModel->has_B_integrand()) {
for (field = fields.begin(); field != fields.end(); field++) {
thisFieldName = field->first;
int index = (int) thisFieldName;
if ( rhs_mask(index,FLUX) ) {
B_fluxes[thisFieldName].reserve(nsd);
}
}
}
// evaluate physics model
if (singleMaterial)
{
if(physicsModel->has_N_integrand()) {
physicsModel->N_integrand(rhs_mask, fieldsAtIPs, grad_fieldsAtIPs,
N_fluxes);
}
if(physicsModel->has_B_integrand()) {
physicsModel->B_integrand(rhs_mask, fieldsAtIPs, grad_fieldsAtIPs,
B_fluxes);
}
}
else
{
// NOTE this copy in fields/copy out fluxes is extremely slow
// NOTE alternates:
// 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
FIELDS groupN_fluxes, fieldsGroup;
GRAD_FIELDS groupB_fluxes, grad_fieldsGroup;
if(physicsModel->has_B_integrand()) {
for (field = fields.begin(); field != fields.end(); field++) {
thisFieldName = field->first;
int index = (int) thisFieldName;
int ndof = (field->second).nCols();
grad_fieldsGroup[thisFieldName].assign(nsd,DENS_MAT(nips,ndof));
N_fluxes[thisFieldName].reset(nips,ndof);
B_fluxes[thisFieldName].assign(nsd,DENS_MAT(nips,ndof));
if ( rhs_mask(index,FLUX) ) {
groupB_fluxes[thisFieldName].reserve(nsd);
}
}
}
// copy fields
for ( int imat = 0; imat < pointMaterialGroups.get_length(); imat++)
{
const set<int> & indices = pointMaterialGroups(imat);
int npts = indices.size();
int i = 0;
for (set<int>::const_iterator iter=indices.begin();
iter != indices.end(); iter++, i++ ) {
for (field = fields.begin(); field != fields.end(); field++) {
thisFieldName = field->first;
int ndof = (field->second).nCols();
fieldsGroup[thisFieldName].reset(npts,ndof);
for (int j = 0; j < nsd; ++j) {
(grad_fieldsGroup[thisFieldName][j]).reset(npts,ndof);
}
for (int dof = 0; dof < ndof; ++dof) {
fieldsGroup[thisFieldName](i,dof)
= fieldsAtIPs[thisFieldName](*iter,dof);
for (int j = 0; j < nsd; ++j) {
grad_fieldsGroup[thisFieldName][j](i,dof)
= grad_fieldsAtIPs[thisFieldName][j](*iter,dof);
}
}
}
}
// calculate forces, & assemble
if(physicsModel->has_N_integrand()) {
physicsModel->N_integrand(rhs_mask, fieldsGroup, grad_fieldsGroup,
groupN_fluxes, imat);
for (field = fields.begin(); field != fields.end(); field++) {
thisFieldName = field->first;
int index = (int) thisFieldName;
int ndof = (field->second).nCols();
if ( rhs_mask(index,SOURCE) ) {
int i = 0;
for (set<int>::const_iterator iter=indices.begin();
iter != indices.end(); iter++, i++ ) {
for (int dof = 0; dof < ndof; ++dof) {
N_fluxes[thisFieldName](*iter,dof)
= groupN_fluxes[thisFieldName](i,dof);
}
}
}
}
}
// calculate forces, & assemble
if(physicsModel->has_B_integrand()) {
physicsModel->B_integrand(rhs_mask, fieldsGroup, grad_fieldsGroup,
groupB_fluxes, imat);
for (field = fields.begin(); field != fields.end(); field++) {
thisFieldName = field->first;
int index = (int) thisFieldName;
int ndof = (field->second).nCols();
if ( rhs_mask(index,FLUX) ) {
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) {
B_fluxes[thisFieldName][j](*iter,dof)
= groupB_fluxes[thisFieldName][j](i,dof);
}
}
}
}
}
}
}
}
// assemble
for (field = fields.begin(); field != fields.end(); field++)
{
thisFieldName = field->first;
int index = (int) thisFieldName;
if(physicsModel->has_N_integrand()) {
if ( rhs_mask(index,SOURCE) ) {
rhs_local[thisFieldName]
+= N.transMat(weights*N_fluxes[thisFieldName]);
}
}
if(physicsModel->has_B_integrand()) {
if ( rhs_mask(index,FLUX) ) {
for (int j = 0; j < nsd; ++j) {
rhs_local[thisFieldName]
+= dN[j].transMat(weights*B_fluxes[thisFieldName][j]);
}
}
}
}
}
// Share information between processors
for (field = rhs.begin(); field != rhs.end(); field++) {
thisFieldName = field->first;
if (rhs_mask(thisFieldName,FLUX) || rhs_mask(thisFieldName,SOURCE)) {
int count = rhs_local[thisFieldName].size();
//rhs_local[thisFieldName].print("RHS LCL");
LammpsInterface::instance()->allsum(
rhs_local[thisFieldName].get_ptr(),
rhs[thisFieldName].get_ptr(), count);
}
}
}
// -------------------------------------------------------------
// compute flux for post processing
// -------------------------------------------------------------
void FE_Engine::compute_flux(const Array2D<bool> &rhs_mask,
const FIELDS &fields,
const PhysicsModel * physicsModel,
const Array<int> & elementMaterials,
GRAD_FIELDS &flux,
const Array<bool> *element_mask) const
{
//NOTE move to initialization and make workspace
const int nNodesPerElement = feMesh_->get_nNodesPerElement();
const int nIPsPerElement = feMesh_->get_nIPsPerElement();
const int nsd = feMesh_->get_nSpatialDimensions();
const int nelems = feMesh_->get_nElements();
DENS_MAT Nmat;
//NOTE move to initialization and make workspace
DENS_MAT N(nIPsPerElement,nNodesPerElement);
vector<DENS_MAT> dN(nsd);
dN.assign(nsd, DENS_MAT(nIPsPerElement,nNodesPerElement));
FIELDS fieldsAtIPs, localElementFields;
GRAD_FIELDS grad_fieldsAtIPs;
FIELDS::const_iterator field;
FieldName thisFieldName;
int numFieldDOF;
for (field = fields.begin(); field != fields.end(); field++) {
thisFieldName = field->first;
if (rhs_mask(thisFieldName,FLUX)) {
int nrows = (field->second).nRows();
int ncols = (field->second).nCols();
flux[thisFieldName].assign(nsd, DENS_MAT(nrows,ncols));
}
}
DIAG_MAT weights(nIPsPerElement,nIPsPerElement);
Array<int> conn(nNodesPerElement);
// element wise assembly
int inode = -1;
for (int ielem=0; ielem < nelems; ielem++)
{
// if element is masked, skip it
if (element_mask && !(*element_mask)(ielem)) continue;
// material id
int imat = elementMaterials(ielem);
// evaluate shape functions
feMesh_->shape_function(ielem, N, dN, weights);
int nips = weights.nRows();
// get connectivity
feMesh_->element_connectivity_unique(ielem, conn);
// compute fields and gradients of fields ips of this element
for (field = fields.begin(); field != fields.end(); field++)
{
// field values at all nodes
const DENS_MAT &vI = field->second;
// field values at element nodes
DENS_MAT &vIe = localElementFields[field->first];
// field values at integration points -> to be computed
DENS_MAT &vP = fieldsAtIPs[field->first];
// gradients of field at integration points -> to be computed
vector<DENS_MAT> &dvP = grad_fieldsAtIPs[field->first];
numFieldDOF = vI.nCols();
vIe.reset(nNodesPerElement, numFieldDOF);
// gather local field
for (int i = 0; i < nNodesPerElement; i++)
for (int j = 0; j < numFieldDOF; j++)
vIe(i,j) = vI(conn(i),j);
// interpolate field at integration points
vP = N*vIe;
dvP.assign(nsd, DENS_MAT(nIPsPerElement, numFieldDOF));
for (int j = 0; j < nsd; ++j) dvP[j] = dN[j]*vIe;
}
// NOTE move scaling of N by weights up here
// evaluate Physics model
// B_fluxes is a set of field gradients
if (physicsModel->has_B_integrand()) {
GRAD_FIELDS B_fluxes;
physicsModel->B_integrand(rhs_mask, fieldsAtIPs, grad_fieldsAtIPs, B_fluxes, imat);
// assemble
for (field = fields.begin(); field != fields.end(); field++)
{
thisFieldName = field->first;
numFieldDOF = (field->second).nCols();
if (!rhs_mask(thisFieldName,FLUX)) continue;
for (int j = 0; j < nsd; j++)
{
Nmat = N.transMat(weights*B_fluxes[thisFieldName][j]);
for (int i = 0; i < nNodesPerElement; ++i)
{
inode = conn(i);
for (int k = 0; k < numFieldDOF; ++k)
{
flux[thisFieldName][j](inode,k) += Nmat(i,k);
} // scatter k
} // scatter i
} // loop over nsd
} //loop on fields
} // if B_integrand
} // element loop
} // function call
//-----------------------------------------------------------------
// boundary flux using native quadrature
//-----------------------------------------------------------------
void FE_Engine::compute_boundary_flux(
const Array2D<bool> & rhs_mask,
const FIELDS & fields,
const PhysicsModel * physicsModel,
const Array<int> & elementMaterials,
const set< pair <int,int> > & faceSet,
FIELDS & rhs) const
{
//NOTE move to initialization and make workspace
int nNodesPerElement = feMesh_->get_nNodesPerElement();
int nIPsPerFace = feMesh_->get_nIPsPerFace();
int nsd = feMesh_->get_nSpatialDimensions();
int nelems = feMesh_->get_nElements();
DENS_MAT Nmat(nNodesPerElement,3);
FieldName thisFieldName;
FIELDS::const_iterator field;
for (field = fields.begin(); field != fields.end(); field++) {
thisFieldName = field->first;
if (rhs_mask(thisFieldName,FLUX)) {
int nrows = (field->second).nRows();
int ncols = (field->second).nCols();
rhs [thisFieldName].reset(nrows,ncols);
}
}
//NOTE move to initialization and make workspace
// sizing working arrays
DENS_MAT N(nIPsPerFace,nNodesPerElement);
vector< DENS_MAT > dN(nsd), nN(nsd);
for (int i = 0; i < nsd; i++) {
dN[i].reset(nIPsPerFace,nNodesPerElement);
nN[i].reset(nIPsPerFace,nNodesPerElement);
}
GRAD_FIELDS grad_fieldsAtIPs;
FIELDS fieldsAtIPs;
FIELDS localElementFields;
int numFieldDOF;
for (field = fields.begin(); field != fields.end(); field++) {
thisFieldName = field->first;
numFieldDOF = field->second.nCols();
grad_fieldsAtIPs[thisFieldName].reserve(nsd);
for (int i = 0; i < nsd; ++i) {
grad_fieldsAtIPs[thisFieldName].push_back(DENS_MAT(nIPsPerFace,numFieldDOF));
}
fieldsAtIPs[thisFieldName].reset(nIPsPerFace,numFieldDOF);
localElementFields[thisFieldName].reset(nNodesPerElement,numFieldDOF);
}
DiagonalMatrix<double> weights(nIPsPerFace,nIPsPerFace);
Array<int> conn(nNodesPerElement);
// element wise assembly
int inode = -1;
set< pair <int,int> >::iterator iter;
for (iter = faceSet.begin(); iter != faceSet.end(); iter++) {
// evaluate shape functions at ips
feMesh_->face_shape_function(*iter, N, dN, nN, weights);
int nips = weights.nRows();
// get connectivity
int ielem = iter->first;
int imat = elementMaterials(ielem);
feMesh_->element_connectivity_unique(ielem, conn);
// interpolate fields and gradients of fields ips of this element
for (field = fields.begin(); field != fields.end(); field++) {
thisFieldName = field->first;
const DenseMatrix<double> * thisField
= (const DENS_MAT *) &(field->second);
numFieldDOF = thisField->nCols();
//NOTE an hopefully unnecessary copy (should alias)
for (int i = 0; i < nNodesPerElement; i++) {
for (int j = 0; j < numFieldDOF; j++) {
localElementFields[thisFieldName](i,j) = (*thisField)(conn(i),j);
}
}
// ips X dof = ips X nodes * nodes X dof
fieldsAtIPs[thisFieldName] = N*localElementFields[thisFieldName];
for (int j = 0; j < nsd; ++j) {
grad_fieldsAtIPs[thisFieldName][j] = dN[j]*localElementFields[thisFieldName];
}
}
// Evaluate- physics model
// do nothing for N_integrand
// nN : precomputed and held by ATC_Transfer
if(physicsModel->has_B_integrand()) {
map<FieldName, vector< DENS_MAT > > B_fluxes;
for (field = fields.begin(); field != fields.end(); field++) {
thisFieldName = field->first;
if ( rhs_mask(thisFieldName,FLUX) ) {
B_fluxes[thisFieldName].reserve(nsd);
}
}
physicsModel->B_integrand(rhs_mask, fieldsAtIPs, grad_fieldsAtIPs, B_fluxes,imat);
// assemble
for (field = fields.begin(); field != fields.end(); field++) {
thisFieldName = field->first;
int index = (int) thisFieldName;
const DenseMatrix<double> * thisField
= (const DENS_MAT *) &(field->second);
if ( rhs_mask(index,FLUX) ) {
numFieldDOF = thisField->nCols();
Nmat.reset(nNodesPerElement,numFieldDOF);
for (int j = 0; j < nsd; j++) {
// nodes X dof = nodes X ips * ips X dof ???
Nmat = nN[j].transMat(weights*B_fluxes[thisFieldName][j]);
for (int i = 0; i < nNodesPerElement; ++i) {
inode = conn(i);
for (int k = 0; k < numFieldDOF; ++k) {
rhs[thisFieldName](inode,k) += Nmat(i,k);
}
}
}
}
}
}
}
}
// -------------------------------------------------------------
// compute boundary flux using given quadrature and interpolation
// -------------------------------------------------------------
void FE_Engine::compute_boundary_flux(
const Array2D<bool> & rhs_mask,
const FIELDS & fields,
const PhysicsModel * physicsModel,
const Array< int > & elementMaterials,
const Array< set<int> > & pointMaterialGroups,
const DIAG_MAT & weights,
const SPAR_MAT & N,
const GRAD_SHPFCN & dN,
const DIAG_MAT & flux_mask,
FIELDS & flux) const
{
int nNodesUnique = feMesh_->get_nNodesUnique();
FieldName thisFieldName;
map<FieldName, DENS_MAT > rhs;
map<FieldName, DENS_MAT > rhsSubset;
map<FieldName, DENS_MAT > rhsSubset_local;
FIELDS::const_iterator field;
for (field = fields.begin(); field != fields.end(); field++) {
thisFieldName = field->first;
if (rhs_mask(thisFieldName,FLUX)) {
int nrows = (field->second).nRows();
int ncols = (field->second).nCols();
rhs [thisFieldName].reset(nrows,ncols);
rhsSubset[thisFieldName].reset(nrows,ncols);
}
}
// R_I = - int_Omega Delta N_I .q dV
compute_rhs_vector(rhs_mask, fields, physicsModel, elementMaterials, rhs);
// R_I^md = - int_Omega^md Delta N_I .q dV
compute_rhs_vector(rhs_mask, 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 < NUM_FIELDS; ++n) {
FieldName thisFieldName = (FieldName) n;
if ( rhs_mask(thisFieldName,FLUX) ) {
flux[thisFieldName]
= rhsSubset[thisFieldName] - flux_mask*rhs[thisFieldName];
}
}
}
//-----------------------------------------------------------------
// 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
{
int nNodesPerElement = feMesh_->get_nNodesPerElement();
int nNodesPerFace = feMesh_->get_nNodesPerFace();
int nIPsPerFace = feMesh_->get_nIPsPerFace();
int nsd = feMesh_->get_nSpatialDimensions();
FieldName thisField;
int nFieldDOF = 1;
XT_Function * f = NULL;
//NOTE move arrays to initialization and make workspace
// sizing working arrays
vector<DENS_MAT> dN, nN;
dN.assign(nsd, DENS_MAT(nIPsPerFace, nNodesPerElement));
nN.assign(nsd, DENS_MAT(nIPsPerFace, nNodesPerElement));
DENS_MAT N(nIPsPerFace,nNodesPerElement);
DENS_MAT xCoords(nsd,nNodesPerElement);
DENS_MAT xAtIPs(nsd,nIPsPerFace);
DENS_MAT Nmat(nNodesPerElement,nsd);
DIAG_MAT weights(nIPsPerFace,nIPsPerFace);
Array<int> conn(nNodesPerElement);
DENS_MAT faceSource(nIPsPerFace,nFieldDOF);
// element wise assembly
SURFACE_SOURCE::const_iterator src_iter;
for (src_iter=sourceFunctions.begin(); src_iter!=sourceFunctions.end(); src_iter++)
{
FieldName thisFieldName = src_iter->first;
if (!fieldMask((int)thisFieldName)) 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;
const Array <XT_Function*> &fs = fset_iter->second;
feMesh_->element_connectivity_unique(elem, conn);
// evaluate location at ips
feMesh_->face_shape_function(face, N, dN, nN, weights);
feMesh_->element_coordinates(elem, xCoords);
MultAB(xCoords,N,xAtIPs,0,1); //xAtIPs = xCoords*(N.transpose());
// interpolate prescribed flux at ips of this element
// NOTE this seems totally redundant
FSET::const_iterator face_iter = fset->find(face);
if (face_iter == fset->end())
{
cout << "face not found" << std::endl;
// NOTE: do something there
}
// NOTE why can't I know the size earlier??
// NOTE only normal flux here
nFieldDOF = (face_iter->second).size();
faceSource.reset(nIPsPerFace,nFieldDOF);
for (int ip = 0; ip < nIPsPerFace; ++ip)
{
for (int idof = 0; idof<nFieldDOF; ++idof)
{
f = fs(idof);
if (!f) continue;
faceSource(ip,idof) = f->f(column(xAtIPs,ip).get_ptr(),time);
}
}
// assemble
Nmat = N.transMat(weights*faceSource);
for (int i = 0; i < nNodesPerElement; ++i)
{
int inode = conn(i);
for (int idof = 0; idof < nFieldDOF; ++idof)
{
nodalSources[thisFieldName](inode,idof) += Nmat(i,idof);
} // end assemble nFieldDOF
} // end assemble nNodesPerElement
} // end fset loop
} // field loop
}
//-----------------------------------------------------------------
// 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
{
int nNodes = feMesh_->get_nNodes();
int nNodesPerElement = feMesh_->get_nNodesPerElement();
int nIPsPerElement = feMesh_->get_nIPsPerElement();
int nsd = feMesh_->get_nSpatialDimensions();
int nelems = feMesh_->get_nElements();
int nFieldDOF = 1;
XT_Function * f = NULL;
//NOTE move arrays to initialization and make workspace
// sizing working arrays
DENS_MAT elemSource(nIPsPerElement,nFieldDOF);
DENS_MAT N(nIPsPerElement,nNodesPerElement);
DENS_MAT xCoords(nsd,nNodesPerElement);
DENS_MAT xAtIPs(nsd,nIPsPerElement);
DENS_MAT Nwg(nNodesPerElement,1);
DENS_MAT Nmat(nNodesPerElement,nsd);
DiagonalMatrix<double> weights(nIPsPerElement,nIPsPerElement);
Array<int> conn(nNodesPerElement);
// element wise assembly
for (int ielem = 0; ielem < nelems; ++ielem) {
feMesh_->element_connectivity_unique(ielem, conn);
// evaluate location at ips
feMesh_->shape_function(ielem, N, weights);
feMesh_->element_coordinates(ielem, xCoords);
xAtIPs =xCoords*(N.transpose());
VOLUME_SOURCE ::const_iterator src_iter;
for (src_iter = sourceFunctions.begin();
src_iter != sourceFunctions.end(); src_iter++) {
FieldName thisField = src_iter->first;
int index = (int) thisField;
if ( fieldMask(index) ) {
const Array2D<XT_Function *> * thisSource
= (const Array2D<XT_Function *> *) &(src_iter->second);
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) {
f = (*thisSource)(ielem,idof);
if (f) {
elemSource(ip,idof) = f->f(column(xAtIPs,ip).get_ptr(),time);
}
}
}
// assemble
Nmat.reset(nNodesPerElement,nFieldDOF);
Nmat = N.transMat(weights*elemSource);
for (int i = 0; i < nNodesPerElement; ++i) {
int inode = conn(i);
for (int idof = 0; idof < nFieldDOF; ++idof) {
nodalSources[thisField](inode,idof) += Nmat(i,idof);
}
}
}
}
}
}
//-----------------------------------------------------------------
// 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 nNodesPerElement = feMesh_->get_nNodesPerElement();
int nIPsPerFace = feMesh_->get_nIPsPerFace();
int nsd = feMesh_->get_nSpatialDimensions();
int nelems = feMesh_->get_nElements();
int numFieldDOF = field.nCols();
double a[3] = {0,0,0};
a[axis] = 1;
// sizing working arrays
DENS_MAT N(nIPsPerFace,nNodesPerElement);
DENS_MAT n(nsd,nIPsPerFace);
DENS_MAT fieldsAtIPs(nIPsPerFace,numFieldDOF);
DENS_MAT localElementFields(nNodesPerElement,numFieldDOF);
DiagonalMatrix<double> weights(nIPsPerFace,nIPsPerFace);
Array<int> conn(nNodesPerElement);
DENS_MAT Nmat(nsd,numFieldDOF);
DENS_MAT integrals(numFieldDOF,nsd);
// element wise assembly
set< pair <int,int> >::iterator iter;
for (iter = faceSet.begin(); iter != faceSet.end(); iter++) {
// evaluate shape functions at ips
feMesh_->face_shape_function(*iter, N, n, weights);
// 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];
}
}
int nips = weights.nRows();
// get connectivity
int ielem = iter->first;
feMesh_->element_connectivity_unique(ielem, conn);
// interpolate fields and gradients of fields ips of this element
for (int i = 0; i < nNodesPerElement; i++) {
for (int j = 0; j < numFieldDOF; 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*weights*fieldsAtIPs;
for (int j = 0; j < nsd; j++) {
for (int k = 0; k < numFieldDOF; ++k) {
integrals(k,j) += Nmat(j,k);
}
}
}
// (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 (numFieldDOF == 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 (numFieldDOF == 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 (numFieldDOF == 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 (numFieldDOF == 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 += tostring(numFieldDOF);
throw ATC_Error(0,msg);
}
}
//-----------------------------------------------------------------
// evaluate shape functions at given points
//-----------------------------------------------------------------
void FE_Engine::evaluate_shape_functions(
const MATRIX & pt_coords,
SPAR_MAT & N,
Array<int> & pointToEltMap) const
{
// Get shape function and derivatives at atomic locations
int npe = feMesh_->get_nNodesPerElement();
int nsd = feMesh_->get_nSpatialDimensions();
int nnodes = feMesh_->get_nNodesUnique();
int npts = pt_coords.nCols();
pointToEltMap.reset(npts);
// loop over point (atom) coordinates
DENS_VEC x(nsd);
Array<int> node_index(npe);
DENS_VEC shp(npe);
N.reset(npts,nnodes);
for (int i = 0; i < npts; ++i) {
for (int k = 0; k < nsd; ++k) {
x(k) = pt_coords(k,i);
}
int eltID;
feMesh_->shape_functions(x,shp,eltID,node_index);
for (int j = 0; j < npe; ++j) {
int jnode = node_index(j);
N.add(i,jnode,shp(j));
}
pointToEltMap(i) = eltID;
}
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,
GRAD_SHPFCN & dN,
Array<int> & pointToEltMap) const
{
// Get shape function and derivatives at atomic locations
int npe = feMesh_->get_nNodesPerElement();
int nsd = feMesh_->get_nSpatialDimensions();
int nnodes = feMesh_->get_nNodesUnique();
int npts = pt_coords.nCols();
pointToEltMap.reset(npts);
// loop over point (atom) coordinates
DENS_VEC x(nsd);
Array<int> node_index(npe);
DENS_VEC shp(npe);
DENS_MAT dshp(nsd,npe);
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(k,i);
}
int eltID;
feMesh_->shape_functions(x,shp,dshp,eltID,node_index);
for (int j = 0; j < npe; ++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));
}
}
pointToEltMap(i) = eltID;
}
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 VECTOR & x,
Array<int>& node_index,
DENS_VEC& shp,
DENS_MAT& dshp,
int & eltID) const
{
// Get shape function and derivatives at a specific point
int nsd = feMesh_->get_nSpatialDimensions();
int npe = feMesh_->get_nNodesPerElement();
dshp.reset(nsd,npe);
feMesh_->shape_functions(x,shp,dshp,eltID,node_index);
}
//-----------------------------------------------------------------
// create a uniform rectangular mesh on a rectangular region
//-----------------------------------------------------------------
void FE_Engine::create_mesh(int nx, int ny, int nz, char * regionName,
int xperiodic, int yperiodic, int zperiodic)
{
double xmin, xmax, ymin, ymax, zmin, zmax;
double xscale, yscale, zscale;
double boxxlo, boxxhi, boxylo, boxyhi, boxzlo, boxzhi;
// check to see if region exists and is it a box, and if so get the bounds
bool isBox;
isBox = atcTransfer_->get_region_bounds(regionName,
xmin, xmax,
ymin, ymax,
zmin, zmax,
xscale,
yscale,
zscale);
if (!isBox) throw ATC_Error(0,"Region for FE mesh is not a box");
feMesh_ = new FE_Uniform3DMesh(nx,ny,nz,
xmin,xmax,
ymin,ymax,
zmin,zmax,
xscale,
yscale,
zscale,
xperiodic,
yperiodic,
zperiodic);
if (LammpsInterface::instance()->comm_rank()==0) {
printf(" ATC:: created FEM Mesh with %i Global Nodes, %i Unique Nodes, and %i Elements\n",
feMesh_->get_nNodes(),feMesh_->get_nNodesUnique(),feMesh_->get_nElements()); }
}
};
c4science · Help