v0.14.0
HookeElement.cpp

Operators and data structures for linear elastic analysisSee as well header file HookeElement.hpp

Implemention of operators for Hooke material. Implementation is extended to the case when the mesh is moving as results of topological changes, also the calculation of material forces and associated tangent matrices are added to implementation.

In other words spatial deformation is small but topological changes large.

/** \file HookeElement.cpp
* \example HookeElement.cpp
* \brief Operators and data structures for linear elastic analysis
*
* See as well header file HookeElement.hpp
*
* Implemention of operators for Hooke material. Implementation is extended to
* the case when the mesh is moving as results of topological changes, also the
* calculation of material forces and associated tangent matrices are added to
* implementation.
*
* In other words spatial deformation is small but topological changes large.
*/
#include <MoFEM.hpp>
using namespace MoFEM;
#include <HookeElement.hpp>
HookeElement::OpCalculateStrainAle::OpCalculateStrainAle(
const std::string row_field, const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
: VolUserDataOperator(row_field, col_field, OPROW, false),
dataAtPts(data_at_pts) {
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
MoFEMErrorCode HookeElement::OpCalculateStrainAle::doWork(int row_side,
EntityType row_type,
EntData &row_data) {
// get number of integration points
const int nb_integration_pts = getGaussPts().size2();
auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
auto t_H = getFTensor2FromMat<3, 3>(*dataAtPts->HMat);
dataAtPts->detHVec->resize(nb_integration_pts, false);
dataAtPts->invHMat->resize(9, nb_integration_pts, false);
dataAtPts->FMat->resize(9, nb_integration_pts, false);
dataAtPts->smallStrainMat->resize(6, nb_integration_pts, false);
auto t_detH = getFTensor0FromVec(*dataAtPts->detHVec);
auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
auto t_F = getFTensor2FromMat<3, 3>(*dataAtPts->FMat);
auto t_strain = getFTensor2SymmetricFromMat<3>(*dataAtPts->smallStrainMat);
for (int gg = 0; gg != nb_integration_pts; ++gg) {
CHKERR invertTensor3by3(t_H, t_detH, t_invH);
t_F(i, j) = t_h(i, k) * t_invH(k, j);
t_strain(i, j) = (t_F(i, j) || t_F(j, i)) / 2.;
t_strain(0, 0) -= 1;
t_strain(1, 1) -= 1;
t_strain(2, 2) -= 1;
++t_strain;
++t_h;
++t_H;
++t_detH;
++t_invH;
++t_F;
}
}
HookeElement::OpCalculateEnergy::OpCalculateEnergy(
const std::string row_field, const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
SmartPetscObj<Vec> ghost_vec)
: VolUserDataOperator(row_field, col_field, OPROW, true),
dataAtPts(data_at_pts), ghostVec(ghost_vec, true) {
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
MoFEMErrorCode HookeElement::OpCalculateEnergy::doWork(int row_side,
EntityType row_type,
EntData &row_data) {
// get number of integration points
const int nb_integration_pts = getGaussPts().size2();
auto t_strain = getFTensor2SymmetricFromMat<3>(*(dataAtPts->smallStrainMat));
auto t_cauchy_stress =
getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
dataAtPts->energyVec->resize(nb_integration_pts, false);
&*(dataAtPts->energyVec->data().begin()));
for (int gg = 0; gg != nb_integration_pts; ++gg) {
t_energy = (t_strain(i, j) * t_cauchy_stress(i, j)) / 2.;
++t_strain;
++t_cauchy_stress;
++t_energy;
}
if (ghostVec.get()) {
// get element volume
double vol = getVolume();
// get intergrayion weights
auto t_w = getFTensor0IntegrationWeight();
auto &det_H = *dataAtPts->detHVec;
&*(dataAtPts->energyVec->data().begin()));
double energy = 0;
for (int gg = 0; gg != nb_integration_pts; ++gg) {
// calculate scalar weight times element volume
double a = t_w * vol;
if (det_H.size()) {
a *= det_H[gg];
}
energy += a * t_energy;
++t_energy;
++t_w;
}
CHKERR VecSetValue(ghostVec, 0, energy, ADD_VALUES);
}
}
HookeElement::OpCalculateEshelbyStress::OpCalculateEshelbyStress(
const std::string row_field, const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> data_at_pts)
: VolUserDataOperator(row_field, col_field, OPROW, true),
dataAtPts(data_at_pts) {
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
MoFEMErrorCode HookeElement::OpCalculateEshelbyStress::doWork(
int row_side, EntityType row_type, EntData &row_data) {
// get number of integration points
const int nb_integration_pts = getGaussPts().size2();
&*(dataAtPts->energyVec->data().begin()));
auto t_cauchy_stress =
getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
auto t_F = getFTensor2FromMat<3, 3>(*(dataAtPts->FMat));
dataAtPts->eshelbyStressMat->resize(9, nb_integration_pts, false);
auto t_eshelby_stress =
getFTensor2FromMat<3, 3>(*(dataAtPts->eshelbyStressMat));
for (int gg = 0; gg != nb_integration_pts; ++gg) {
t_eshelby_stress(i, j) = -t_F(k, i) * t_cauchy_stress(k, j);
t_eshelby_stress(0, 0) += t_energy;
t_eshelby_stress(1, 1) += t_energy;
t_eshelby_stress(2, 2) += t_energy;
++t_cauchy_stress;
++t_energy;
++t_eshelby_stress;
++t_F;
}
}
HookeElement::OpAssemble::OpAssemble(
const std::string row_field, const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> data_at_pts, const char type,
bool symm)
: VolUserDataOperator(row_field, col_field, type, symm),
dataAtPts(data_at_pts) {}
MoFEMErrorCode HookeElement::OpAssemble::doWork(int row_side, int col_side,
EntityType row_type,
EntityType col_type,
EntData &row_data,
EntData &col_data) {
// get number of dofs on row
nbRows = row_data.getIndices().size();
// if no dofs on row, exit that work, nothing to do here
if (!nbRows)
// get number of dofs on column
nbCols = col_data.getIndices().size();
// if no dofs on Columbia, exit nothing to do here
if (!nbCols)
// K_ij matrix will have 3 times the number of degrees of freedom of the
// i-th entity set (nbRows)
// and 3 times the number of degrees of freedom of the j-th entity set
// (nbCols)
K.resize(nbRows, nbCols, false);
K.clear();
// get number of integration points
nbIntegrationPts = getGaussPts().size2();
// check if entity block is on matrix diagonal
if (row_side == col_side && row_type == col_type) {
isDiag = true;
} else {
isDiag = false;
}
// integrate local matrix for entity block
CHKERR iNtegrate(row_data, col_data);
// assemble local matrix
CHKERR aSsemble(row_data, col_data);
}
MoFEMErrorCode HookeElement::OpAssemble::doWork(int row_side,
EntityType row_type,
EntData &row_data) {
// get number of dofs on row
nbRows = row_data.getIndices().size();
// if no dofs on row, exit that work, nothing to do here
if (!nbRows)
nF.resize(nbRows, false);
nF.clear();
// get number of integration points
nbIntegrationPts = getGaussPts().size2();
// integrate local matrix for entity block
CHKERR iNtegrate(row_data);
// assemble local matrix
CHKERR aSsemble(row_data);
}
MoFEMErrorCode HookeElement::OpAssemble::iNtegrate(EntData &row_data,
EntData &col_data) {
SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
};
MoFEMErrorCode HookeElement::OpAssemble::iNtegrate(EntData &row_data) {
SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
};
MoFEMErrorCode HookeElement::OpAssemble::aSsemble(EntData &row_data,
EntData &col_data) {
// get pointer to first global index on row
const int *row_indices = &*row_data.getIndices().data().begin();
// get pointer to first global index on column
const int *col_indices = &*col_data.getIndices().data().begin();
auto &data = *dataAtPts;
if (!data.forcesOnlyOnEntitiesRow.empty()) {
rowIndices.resize(nbRows, false);
noalias(rowIndices) = row_data.getIndices();
row_indices = &rowIndices[0];
VectorDofs &dofs = row_data.getFieldDofs();
VectorDofs::iterator dit = dofs.begin();
for (int ii = 0; dit != dofs.end(); dit++, ii++) {
if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
data.forcesOnlyOnEntitiesRow.end()) {
rowIndices[ii] = -1;
}
}
}
if (!data.forcesOnlyOnEntitiesCol.empty()) {
colIndices.resize(nbCols, false);
noalias(colIndices) = col_data.getIndices();
col_indices = &colIndices[0];
VectorDofs &dofs = col_data.getFieldDofs();
VectorDofs::iterator dit = dofs.begin();
for (int ii = 0; dit != dofs.end(); dit++, ii++) {
if (data.forcesOnlyOnEntitiesCol.find((*dit)->getEnt()) ==
data.forcesOnlyOnEntitiesCol.end()) {
colIndices[ii] = -1;
}
}
}
Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
: getFEMethod()->snes_B;
// assemble local matrix
CHKERR MatSetValues(B, nbRows, row_indices, nbCols, col_indices,
&*K.data().begin(), ADD_VALUES);
if (!isDiag && sYmm) {
// if not diagonal term and since global matrix is symmetric assemble
// transpose term.
transK.resize(K.size2(), K.size1(), false);
noalias(transK) = trans(K);
CHKERR MatSetValues(B, nbCols, col_indices, nbRows, row_indices,
&*transK.data().begin(), ADD_VALUES);
}
}
MoFEMErrorCode HookeElement::OpAssemble::aSsemble(EntData &row_data) {
// get pointer to first global index on row
const int *row_indices = &*row_data.getIndices().data().begin();
auto &data = *dataAtPts;
if (!data.forcesOnlyOnEntitiesRow.empty()) {
rowIndices.resize(nbRows, false);
noalias(rowIndices) = row_data.getIndices();
row_indices = &rowIndices[0];
VectorDofs &dofs = row_data.getFieldDofs();
VectorDofs::iterator dit = dofs.begin();
for (int ii = 0; dit != dofs.end(); dit++, ii++) {
if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
data.forcesOnlyOnEntitiesRow.end()) {
rowIndices[ii] = -1;
}
}
}
Vec F = getFEMethod()->snes_f;
// assemble local matrix
CHKERR VecSetValues(F, nbRows, row_indices, &*nF.data().begin(), ADD_VALUES);
}
HookeElement::OpRhs_dx::OpRhs_dx(
const std::string row_field, const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
: OpAssemble(row_field, col_field, data_at_pts, OPROW) {}
MoFEMErrorCode HookeElement::OpRhs_dx::iNtegrate(EntData &row_data) {
auto get_tensor1 = [](VectorDouble &v, const int r) {
&v(r + 0), &v(r + 1), &v(r + 2));
};
// get element volume
double vol = getVolume();
// get intergrayion weights
auto t_w = getFTensor0IntegrationWeight();
// get derivatives of base functions on rows
auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
const int row_nb_base_fun = row_data.getN().size2();
auto t_cauchy_stress =
getFTensor2SymmetricFromMat<3>(*dataAtPts->cauchyStressMat);
// iterate over integration points
for (int gg = 0; gg != nbIntegrationPts; ++gg) {
// calculate scalar weight times element volume
double a = t_w * vol;
auto t_nf = get_tensor1(nF, 0);
int rr = 0;
for (; rr != nbRows / 3; ++rr) {
t_nf(i) += a * t_row_diff_base(j) * t_cauchy_stress(i, j);
++t_row_diff_base;
++t_nf;
}
for (; rr != row_nb_base_fun; ++rr)
++t_row_diff_base;
++t_w;
++t_cauchy_stress;
}
}
HookeElement::OpAleRhs_dx::OpAleRhs_dx(
const std::string row_field, const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
: OpAssemble(row_field, col_field, data_at_pts, OPROW) {}
MoFEMErrorCode HookeElement::OpAleRhs_dx::iNtegrate(EntData &row_data) {
auto get_tensor1 = [](VectorDouble &v, const int r) {
&v(r + 0), &v(r + 1), &v(r + 2));
};
// get element volume
double vol = getVolume();
// get intergrayion weights
auto t_w = getFTensor0IntegrationWeight();
// get derivatives of base functions on rows
auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
const int row_nb_base_fun = row_data.getN().size2();
auto t_cauchy_stress =
getFTensor2SymmetricFromMat<3>(*dataAtPts->cauchyStressMat);
auto &det_H = *dataAtPts->detHVec;
auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
// iterate over integration points
for (int gg = 0; gg != nbIntegrationPts; ++gg) {
// calculate scalar weight times element volume
double a = t_w * vol * det_H[gg];
auto t_nf = get_tensor1(nF, 0);
int rr = 0;
for (; rr != nbRows / 3; ++rr) {
FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
t_nf(i) += a * t_row_diff_base_pulled(j) * t_cauchy_stress(i, j);
++t_row_diff_base;
++t_nf;
}
for (; rr != row_nb_base_fun; ++rr)
++t_row_diff_base;
++t_w;
++t_cauchy_stress;
++t_invH;
}
}
HookeElement::OpAleRhs_dX::OpAleRhs_dX(
const std::string row_field, const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
: OpAssemble(row_field, col_field, data_at_pts, OPROW) {}
MoFEMErrorCode HookeElement::OpAleRhs_dX::iNtegrate(EntData &row_data) {
auto get_tensor1 = [](VectorDouble &v, const int r) {
&v(r + 0), &v(r + 1), &v(r + 2));
};
// get element volume
double vol = getVolume();
// get intergrayion weights
auto t_w = getFTensor0IntegrationWeight();
// get derivatives of base functions on rows
auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
const int row_nb_base_fun = row_data.getN().size2();
auto t_eshelby_stress =
getFTensor2FromMat<3, 3>(*dataAtPts->eshelbyStressMat);
auto &det_H = *dataAtPts->detHVec;
auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
// iterate over integration points
for (int gg = 0; gg != nbIntegrationPts; ++gg) {
// calculate scalar weight times element volume
double a = t_w * vol * det_H[gg];
auto t_nf = get_tensor1(nF, 0);
int rr = 0;
for (; rr != nbRows / 3; ++rr) {
FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
t_nf(i) += a * t_row_diff_base_pulled(j) * t_eshelby_stress(i, j);
++t_row_diff_base;
++t_nf;
}
for (; rr != row_nb_base_fun; ++rr)
++t_row_diff_base;
++t_w;
++t_eshelby_stress;
++t_invH;
}
}
MoFEM::Interface &m_field,
boost::shared_ptr<map<int, BlockData>> &block_sets_ptr) {
if (!block_sets_ptr)
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Pointer to block of sets is null");
m_field, BLOCKSET | MAT_ELASTICSET, it)) {
Mat_Elastic mydata;
CHKERR it->getAttributeDataStructure(mydata);
int id = it->getMeshsetId();
auto &block_data = (*block_sets_ptr)[id];
EntityHandle meshset = it->getMeshset();
CHKERR m_field.get_moab().get_entities_by_dimension(meshset, 3,
block_data.tEts, true);
block_data.iD = id;
block_data.E = mydata.data.Young;
block_data.PoissonRatio = mydata.data.Poisson;
}
}
MoFEM::Interface &m_field,
boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
const std::string element_name, const std::string x_field,
const std::string X_field, const bool ale) {
if (!block_sets_ptr)
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Pointer to block of sets is null");
CHKERR m_field.add_finite_element(element_name, MF_ZERO);
CHKERR m_field.modify_finite_element_add_field_row(element_name, x_field);
CHKERR m_field.modify_finite_element_add_field_col(element_name, x_field);
CHKERR m_field.modify_finite_element_add_field_data(element_name, x_field);
if (m_field.check_field(X_field)) {
if (ale) {
CHKERR m_field.modify_finite_element_add_field_row(element_name, X_field);
CHKERR m_field.modify_finite_element_add_field_col(element_name, X_field);
}
CHKERR m_field.modify_finite_element_add_field_data(element_name, X_field);
}
for (auto &m : (*block_sets_ptr)) {
CHKERR m_field.add_ents_to_finite_element_by_dim(m.second.tEts, 3,
element_name);
}
}
boost::shared_ptr<ForcesAndSourcesCore> fe_lhs_ptr,
boost::shared_ptr<ForcesAndSourcesCore> fe_rhs_ptr,
boost::shared_ptr<map<int, BlockData>> block_sets_ptr,
const std::string x_field, const std::string X_field, const bool ale,
const bool field_disp, const EntityType type,
boost::shared_ptr<DataAtIntegrationPts> data_at_pts) {
if (!block_sets_ptr)
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Pointer to block of sets is null");
if (!data_at_pts)
data_at_pts = boost::make_shared<DataAtIntegrationPts>();
if (fe_lhs_ptr) {
if (ale == PETSC_FALSE) {
if (type == MBPRISM) {
boost::shared_ptr<MatrixDouble> inv_jac_ptr(new MatrixDouble);
fe_lhs_ptr->getOpPtrVector().push_back(
new OpCalculateInvJacForFatPrism(inv_jac_ptr));
fe_lhs_ptr->getOpPtrVector().push_back(
new OpSetInvJacH1ForFatPrism(inv_jac_ptr));
}
fe_lhs_ptr->getOpPtrVector().push_back(
x_field, x_field, block_sets_ptr, data_at_pts));
fe_lhs_ptr->getOpPtrVector().push_back(
new OpLhs_dx_dx<0>(x_field, x_field, data_at_pts));
} else {
if (type == MBPRISM) {
boost::shared_ptr<MatrixDouble> inv_jac_ptr(new MatrixDouble);
fe_lhs_ptr->getOpPtrVector().push_back(
new OpCalculateInvJacForFatPrism(inv_jac_ptr));
fe_lhs_ptr->getOpPtrVector().push_back(
new OpSetInvJacH1ForFatPrism(inv_jac_ptr));
}
fe_lhs_ptr->getOpPtrVector().push_back(
new OpCalculateVectorFieldGradient<3, 3>(X_field, data_at_pts->HMat));
fe_lhs_ptr->getOpPtrVector().push_back(
x_field, x_field, block_sets_ptr, data_at_pts));
fe_lhs_ptr->getOpPtrVector().push_back(
new OpCalculateVectorFieldGradient<3, 3>(x_field, data_at_pts->hMat));
fe_lhs_ptr->getOpPtrVector().push_back(
new OpCalculateStrainAle(x_field, x_field, data_at_pts));
fe_lhs_ptr->getOpPtrVector().push_back(
new OpCalculateStress<0>(x_field, x_field, data_at_pts));
fe_lhs_ptr->getOpPtrVector().push_back(
new OpAleLhs_dx_dx<0>(x_field, x_field, data_at_pts));
fe_lhs_ptr->getOpPtrVector().push_back(
new OpAleLhs_dx_dX<0>(x_field, X_field, data_at_pts));
fe_lhs_ptr->getOpPtrVector().push_back(
new OpCalculateEnergy(X_field, X_field, data_at_pts));
fe_lhs_ptr->getOpPtrVector().push_back(
new OpCalculateEshelbyStress(X_field, X_field, data_at_pts));
fe_lhs_ptr->getOpPtrVector().push_back(
new OpAleLhs_dX_dX<0>(X_field, X_field, data_at_pts));
fe_lhs_ptr->getOpPtrVector().push_back(
new OpAleLhsPre_dX_dx<0>(X_field, x_field, data_at_pts));
fe_lhs_ptr->getOpPtrVector().push_back(
new OpAleLhs_dX_dx(X_field, x_field, data_at_pts));
}
}
if (fe_rhs_ptr) {
if (ale == PETSC_FALSE) {
if (type == MBPRISM) {
boost::shared_ptr<MatrixDouble> inv_jac_ptr(new MatrixDouble);
fe_rhs_ptr->getOpPtrVector().push_back(
new OpCalculateInvJacForFatPrism(inv_jac_ptr));
fe_rhs_ptr->getOpPtrVector().push_back(
new OpSetInvJacH1ForFatPrism(inv_jac_ptr));
}
fe_rhs_ptr->getOpPtrVector().push_back(
new OpCalculateVectorFieldGradient<3, 3>(x_field, data_at_pts->hMat));
fe_rhs_ptr->getOpPtrVector().push_back(
new OpCalculateHomogeneousStiffness<0>(x_field, x_field,
block_sets_ptr, data_at_pts));
if (field_disp) {
fe_rhs_ptr->getOpPtrVector().push_back(
new OpCalculateStrain<true>(x_field, x_field, data_at_pts));
} else {
fe_rhs_ptr->getOpPtrVector().push_back(
new OpCalculateStrain<false>(x_field, x_field, data_at_pts));
}
fe_rhs_ptr->getOpPtrVector().push_back(
new OpCalculateStress<0>(x_field, x_field, data_at_pts));
fe_rhs_ptr->getOpPtrVector().push_back(
new OpRhs_dx(x_field, x_field, data_at_pts));
} else {
if (type == MBPRISM) {
boost::shared_ptr<MatrixDouble> inv_jac_ptr(new MatrixDouble);
fe_rhs_ptr->getOpPtrVector().push_back(
new OpCalculateInvJacForFatPrism(inv_jac_ptr));
fe_rhs_ptr->getOpPtrVector().push_back(
new OpSetInvJacH1ForFatPrism(inv_jac_ptr));
}
fe_rhs_ptr->getOpPtrVector().push_back(
new OpCalculateVectorFieldGradient<3, 3>(X_field, data_at_pts->HMat));
fe_rhs_ptr->getOpPtrVector().push_back(
new OpCalculateHomogeneousStiffness<0>(x_field, x_field,
block_sets_ptr, data_at_pts));
fe_rhs_ptr->getOpPtrVector().push_back(
new OpCalculateVectorFieldGradient<3, 3>(x_field, data_at_pts->hMat));
fe_rhs_ptr->getOpPtrVector().push_back(
new OpCalculateStrainAle(x_field, x_field, data_at_pts));
fe_rhs_ptr->getOpPtrVector().push_back(
new OpCalculateStress<0>(x_field, x_field, data_at_pts));
fe_rhs_ptr->getOpPtrVector().push_back(
new OpAleRhs_dx(x_field, x_field, data_at_pts));
fe_rhs_ptr->getOpPtrVector().push_back(
new OpCalculateEnergy(X_field, X_field, data_at_pts));
fe_rhs_ptr->getOpPtrVector().push_back(
new OpCalculateEshelbyStress(X_field, X_field, data_at_pts));
fe_rhs_ptr->getOpPtrVector().push_back(
new OpAleRhs_dX(X_field, X_field, data_at_pts));
}
}
}
DM dm, boost::shared_ptr<map<int, BlockData>> block_sets_ptr,
const std::string x_field, const std::string X_field, const bool ale,
const bool field_disp, SmartPetscObj<Vec> &v_energy) {
MoFEM::Interface *m_field_ptr;
CHKERR DMoFEMGetInterfacePtr(dm, &m_field_ptr);
v_energy = createVectorMPI(m_field_ptr->get_comm(), PETSC_DECIDE, 1);
boost::shared_ptr<DataAtIntegrationPts> data_at_pts(
new DataAtIntegrationPts());
auto fe_ptr =
boost::make_shared<VolumeElementForcesAndSourcesCore>(*m_field_ptr);
fe_ptr->getRuleHook = [](const double, const double, const double o) {
return 2 * o;
};
if (m_field_ptr->check_field("MESH_NODE_POSITIONS")) CHKERR addHOOpsVol(
"MESH_NODE_POSITIONS", *fe_ptr, true, false, false, false);
using FatPrismElementForcesAndSourcesCore::
FatPrismElementForcesAndSourcesCore;
int getRuleTrianglesOnly(int order) { return 2 * order; }
int getRuleThroughThickness(int order) { return 2 * order; }
};
boost::shared_ptr<ForcesAndSourcesCore> prism_fe_ptr(
new PrismFE(*m_field_ptr));
auto push_ops = [&](boost::shared_ptr<ForcesAndSourcesCore> fe_ptr,
EntityType type) {
boost::shared_ptr<MatrixDouble> inv_jac_ptr(new MatrixDouble);
if (ale == PETSC_FALSE) {
if (type == MBPRISM) {
fe_ptr->getOpPtrVector().push_back(
new OpCalculateInvJacForFatPrism(inv_jac_ptr));
fe_ptr->getOpPtrVector().push_back(
new OpSetInvJacH1ForFatPrism(inv_jac_ptr));
}
fe_ptr->getOpPtrVector().push_back(
new OpCalculateVectorFieldGradient<3, 3>(x_field, data_at_pts->hMat));
fe_ptr->getOpPtrVector().push_back(new OpCalculateHomogeneousStiffness<0>(
x_field, x_field, block_sets_ptr, data_at_pts));
if (field_disp) {
fe_ptr->getOpPtrVector().push_back(
new OpCalculateStrain<true>(x_field, x_field, data_at_pts));
} else {
fe_ptr->getOpPtrVector().push_back(
new OpCalculateStrain<false>(x_field, x_field, data_at_pts));
}
fe_ptr->getOpPtrVector().push_back(
new OpCalculateStress<0>(x_field, x_field, data_at_pts));
fe_ptr->getOpPtrVector().push_back(
new OpCalculateEnergy(X_field, X_field, data_at_pts, v_energy));
} else {
if (type == MBPRISM) {
fe_ptr->getOpPtrVector().push_back(
new OpCalculateInvJacForFatPrism(inv_jac_ptr));
fe_ptr->getOpPtrVector().push_back(
new OpSetInvJacH1ForFatPrism(inv_jac_ptr));
}
fe_ptr->getOpPtrVector().push_back(
new OpCalculateVectorFieldGradient<3, 3>(X_field, data_at_pts->HMat));
fe_ptr->getOpPtrVector().push_back(new OpCalculateHomogeneousStiffness<0>(
x_field, x_field, block_sets_ptr, data_at_pts));
fe_ptr->getOpPtrVector().push_back(
new OpCalculateVectorFieldGradient<3, 3>(x_field, data_at_pts->hMat));
fe_ptr->getOpPtrVector().push_back(
new OpCalculateStrainAle(x_field, x_field, data_at_pts));
fe_ptr->getOpPtrVector().push_back(
new OpCalculateStress<0>(x_field, x_field, data_at_pts));
fe_ptr->getOpPtrVector().push_back(
new OpCalculateEnergy(X_field, X_field, data_at_pts, v_energy));
}
};
CHKERR push_ops(fe_ptr, MBTET);
CHKERR push_ops(prism_fe_ptr, MBPRISM);
CHKERR VecZeroEntries(v_energy);
fe_ptr->snes_ctx = SnesMethod::CTX_SNESNONE;
CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", fe_ptr);
CHKERR DMoFEMLoopFiniteElements(dm, "ELASTIC", prism_fe_ptr);
CHKERR VecAssemblyBegin(v_energy);
CHKERR VecAssemblyEnd(v_energy);
}
MoFEMErrorCode HookeElement::OpAleLhs_dX_dx::iNtegrate(EntData &row_data,
EntData &col_data) {
// get sub-block (3x3) of local stiffens matrix, here represented by
// second order tensor
auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
&m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
&m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
&m(r + 2, c + 2));
};
// get element volume
double vol = getVolume();
// get intergrayion weights
auto t_w = getFTensor0IntegrationWeight();
// get derivatives of base functions on rows
auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
const int row_nb_base_fun = row_data.getN().size2();
auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
auto &det_H = *dataAtPts->detHVec;
auto get_eshelby_stress_dx = [this]() {
t_eshelby_stress_dx;
int mm = 0;
for (int ii = 0; ii != 3; ++ii)
for (int jj = 0; jj != 3; ++jj)
for (int kk = 0; kk != 3; ++kk)
for (int ll = 0; ll != 3; ++ll)
t_eshelby_stress_dx.ptr(ii, jj, kk, ll) =
&(*dataAtPts->eshelbyStress_dx)(mm++, 0);
return t_eshelby_stress_dx;
};
auto t_eshelby_stress_dx = get_eshelby_stress_dx();
// iterate over integration points
for (int gg = 0; gg != nbIntegrationPts; ++gg) {
// calculate scalar weight times element volume
double a = t_w * vol * det_H[gg];
// iterate over row base functions
int rr = 0;
for (; rr != nbRows / 3; ++rr) {
// get sub matrix for the row
auto t_m = get_tensor2(K, 3 * rr, 0);
FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
t_row_stress_dx(i, k, l) =
a * t_row_diff_base_pulled(j) * t_eshelby_stress_dx(i, j, k, l);
// get derivatives of base functions for columns
auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
// iterate column base functions
for (int cc = 0; cc != nbCols / 3; ++cc) {
t_m(i, k) += t_row_stress_dx(i, k, l) * t_col_diff_base(l);
// move to next column base function
++t_col_diff_base;
// move to next block of local stiffens matrix
++t_m;
}
// move to next row base function
++t_row_diff_base;
}
for (; rr != row_nb_base_fun; ++rr)
++t_row_diff_base;
++t_w;
++t_invH;
++t_eshelby_stress_dx;
}
}
HookeElement::OpAleLhsWithDensity_dX_dX::OpAleLhsWithDensity_dX_dX(
const std::string row_field, const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts, const double rho_n,
const double rho_0)
: OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, false),
rhoAtGaussPtsPtr(rho_at_gauss_pts),
rhoGradAtGaussPtsPtr(rho_grad_at_gauss_pts), rhoN(rho_n), rHo0(rho_0) {}
HookeElement::OpAleLhsWithDensity_dX_dX::iNtegrate(EntData &row_data,
EntData &col_data) {
// get sub-block (3x3) of local stiffens matrix, here represented by
// second order tensor
auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
&m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
&m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
&m(r + 2, c + 2));
};
// get element volume
double vol = getVolume();
// get intergrayion weights
auto t_w = getFTensor0IntegrationWeight();
// get derivatives of base functions on rows
auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
const int row_nb_base_fun = row_data.getN().size2();
// Elastic stiffness tensor (4th rank tensor with minor and major
// symmetry)
auto rho = getFTensor0FromVec(*rhoAtGaussPtsPtr);
auto t_grad_rho = getFTensor1FromMat<3>(*rhoGradAtGaussPtsPtr);
auto t_eshelby_stress =
getFTensor2FromMat<3, 3>(*dataAtPts->eshelbyStressMat);
// auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
auto &det_H = *dataAtPts->detHVec;
// iterate over integration points
for (int gg = 0; gg != nbIntegrationPts; ++gg) {
// calculate scalar weight times element volume
double a = t_w * vol * det_H[gg];
const double stress_dho_coef = (rhoN / rho);
// (rhoN / rHo0) * pow(rho / rHo0, rhoN - 1.) * (1. / pow(rho / rHo0,
// rhoN));
// iterate over row base functions
int rr = 0;
for (; rr != nbRows / 3; ++rr) {
// get sub matrix for the row
auto t_m = get_tensor2(K, 3 * rr, 0);
FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
t_row_stress(i) = a * t_row_diff_base_pulled(j) * t_eshelby_stress(i, j);
// get derivatives of base functions for columns
// auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
auto t_col_base = col_data.getFTensor0N(gg, 0);
// iterate column base functions
for (int cc = 0; cc != nbCols / 3; ++cc) {
t_m(i, k) +=
t_row_stress(i) * stress_dho_coef * t_grad_rho(k) * t_col_base;
// move to next column base function
++t_col_base;
// move to next block of local stiffens matrix
++t_m;
}
// move to next row base function
++t_row_diff_base;
}
for (; rr != row_nb_base_fun; ++rr)
++t_row_diff_base;
// move to next integration weight
++t_w;
++t_eshelby_stress;
++t_invH;
++rho;
++t_grad_rho;
}
}
HookeElement::OpAleLhsWithDensity_dx_dX::OpAleLhsWithDensity_dx_dX(
const std::string row_field, const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts, const double rho_n,
const double rho_0)
: OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, false),
rhoAtGaussPtsPtr(rho_at_gauss_pts),
rhoGradAtGaussPtsPtr(rho_grad_at_gauss_pts), rhoN(rho_n), rHo0(rho_0) {}
HookeElement::OpAleLhsWithDensity_dx_dX::iNtegrate(EntData &row_data,
EntData &col_data) {
// get sub-block (3x3) of local stiffens matrix, here represented by
// second order tensor
auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
&m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
&m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
&m(r + 2, c + 2));
};
// get element volume
double vol = getVolume();
// get integration weights
auto t_w = getFTensor0IntegrationWeight();
// get derivatives of base functions on rows
auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
const int row_nb_base_fun = row_data.getN().size2();
auto rho = getFTensor0FromVec(*rhoAtGaussPtsPtr);
auto t_grad_rho = getFTensor1FromMat<3>(*rhoGradAtGaussPtsPtr);
auto t_cauchy_stress =
getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
// auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
auto &det_H = *dataAtPts->detHVec;
// iterate over integration points
for (int gg = 0; gg != nbIntegrationPts; ++gg) {
// calculate scalar weight times element volume
double a = t_w * vol * det_H[gg];
const double stress_dho_coef = (rhoN / rho);
// (rhoN / rHo0) * pow(rho / rHo0, rhoN - 1.) * (1. / pow(rho / rHo0,
// rhoN)); iterate over row base functions
int rr = 0;
for (; rr != nbRows / 3; ++rr) {
// get sub matrix for the row
auto t_m = get_tensor2(K, 3 * rr, 0);
FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
t_row_stress(i) = a * t_row_diff_base_pulled(j) * t_cauchy_stress(i, j);
// get derivatives of base functions for columns
auto t_col_base = col_data.getFTensor0N(gg, 0);
// iterate column base functions
for (int cc = 0; cc != nbCols / 3; ++cc) {
t_m(i, k) +=
t_row_stress(i) * stress_dho_coef * t_grad_rho(k) * t_col_base;
++t_col_base;
// move to next block of local stiffens matrix
++t_m;
}
// move to next row base function
++t_row_diff_base;
}
for (; rr != row_nb_base_fun; ++rr)
++t_row_diff_base;
// move to next integration weight
++t_w;
// ++t_D;
++t_cauchy_stress;
++t_invH;
// ++t_h;
++rho;
++t_grad_rho;
}
}
HookeElement::OpCalculateStiffnessScaledByDensityField::
OpCalculateStiffnessScaledByDensityField(
const std::string row_field, const std::string col_field,
boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
boost::shared_ptr<VectorDouble> rho_at_gauss_pts, const double rho_n,
const double rho_0)
: VolUserDataOperator(row_field, col_field, OPROW, true),
blockSetsPtr(block_sets_ptr), dataAtPts(data_at_pts),
rhoAtGaussPtsPtr(rho_at_gauss_pts), rhoN(rho_n), rHo0(rho_0) {
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
MoFEMErrorCode HookeElement::OpCalculateStiffnessScaledByDensityField::doWork(
int row_side, EntityType row_type, EntData &row_data) {
if (!rhoAtGaussPtsPtr)
SETERRQ(PETSC_COMM_SELF, 1, "Calculate density with MWLS first.");
for (auto &m : (*blockSetsPtr)) {
if (m.second.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
m.second.tEts.end()) {
continue;
}
const int nb_integration_pts = getGaussPts().size2();
dataAtPts->stiffnessMat->resize(36, nb_integration_pts, false);
MAT_TO_DDG(dataAtPts->stiffnessMat));
const double young = m.second.E;
const double poisson = m.second.PoissonRatio;
auto rho = getFTensor0FromVec(*rhoAtGaussPtsPtr);
// coefficient used in intermediate calculation
const double coefficient = young / ((1 + poisson) * (1 - 2 * poisson));
for (int gg = 0; gg != nb_integration_pts; ++gg) {
t_D(i, j, k, l) = 0.;
t_D(0, 0, 0, 0) = 1 - poisson;
t_D(1, 1, 1, 1) = 1 - poisson;
t_D(2, 2, 2, 2) = 1 - poisson;
t_D(0, 1, 0, 1) = 0.5 * (1 - 2 * poisson);
t_D(0, 2, 0, 2) = 0.5 * (1 - 2 * poisson);
t_D(1, 2, 1, 2) = 0.5 * (1 - 2 * poisson);
t_D(0, 0, 1, 1) = poisson;
t_D(1, 1, 0, 0) = poisson;
t_D(0, 0, 2, 2) = poisson;
t_D(2, 2, 0, 0) = poisson;
t_D(1, 1, 2, 2) = poisson;
t_D(2, 2, 1, 1) = poisson;
// here the coefficient is modified to take density into account for
// porous materials: E(p) = E * (p / p_0)^n
t_D(i, j, k, l) *= coefficient * pow(rho / rHo0, rhoN);
++t_D;
++rho;
}
}
}
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEM::K
VectorDouble K
Definition: Projection10NodeCoordsOnField.cpp:125
setBlocks
static MoFEMErrorCode setBlocks(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData >> &block_sets_ptr)
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:127
MoFEM::MatSetValues
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
Definition: EntitiesFieldData.hpp:1631
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
EntityHandle
OpLhs_dx_dx
Definition: HookeElement.hpp:410
rho
double rho
Definition: plastic.cpp:140
MoFEM::addHOOpsVol
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
Definition: HODataOperators.hpp:764
OpAleLhsPre_dX_dx
Definition: HookeElement.hpp:536
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::CoreInterface::modify_finite_element_add_field_row
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
MoFEM::Mat_Elastic
Elastic material data structure.
Definition: MaterialBlocks.hpp:139
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
Definition: MeshsetsManager.hpp:49
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MoFEM.hpp
calculateEnergy
static MoFEMErrorCode calculateEnergy(DM dm, boost::shared_ptr< map< int, BlockData >> block_sets_ptr, const std::string x_field, const std::string X_field, const bool ale, const bool field_disp, SmartPetscObj< Vec > &v_energy_ptr)
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1576
sdf.r
int r
Definition: sdf.py:8
FTensor::Tensor2
Definition: Tensor2_value.hpp:16
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::OpCalculateInvJacForFatPrism
Calculate inverse of jacobian for face element.
Definition: UserDataOperators.hpp:3150
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::invertTensor3by3
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
Definition: Templates.hpp:1566
MoFEM::EntitiesFieldData::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition: EntitiesFieldData.hpp:1489
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
FTensor::Tensor3
Definition: Tensor3_value.hpp:12
OpCalculateStrainAle
Definition: HookeElement.hpp:130
MoFEM::CoreInterface::add_finite_element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
OpAleRhs_dx
Definition: HookeElement.hpp:425
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::CoreInterface::modify_finite_element_add_field_col
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
MoFEM::OpSetInvJacH1ForFatPrism
Transform local reference derivatives of shape functions to global derivatives.
Definition: UserDataOperators.hpp:3177
OpCalculateHomogeneousStiffness
Definition: HookeElement.hpp:190
a
constexpr double a
Definition: approx_sphere.cpp:30
double
convert.type
type
Definition: convert.py:64
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
MAT_ELASTICSET
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
Definition: definitions.h:159
OpAleRhs_dX
Definition: HookeElement.hpp:512
OpCalculateEshelbyStress
Definition: HookeElement.hpp:177
OpCalculateEnergy
Definition: HookeElement.hpp:164
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition: EntitiesFieldData.hpp:1201
HookeElement.hpp
MoFEM::VectorDofs
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
Definition: EntitiesFieldData.hpp:23
MoFEM::Mat_Elastic::data
_data_ data
Definition: MaterialBlocks.hpp:155
MoFEM::CoreInterface::check_field
virtual bool check_field(const std::string &name) const =0
check if field is in database
PrismFE
Definition: prism_elements_from_surface.cpp:62
MoFEM::EntitiesFieldData::EntData::getFieldDofs
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
Definition: EntitiesFieldData.hpp:1256
FTensor::Tensor4
Definition: Tensor4_value.hpp:18
MoFEM::CoreInterface::add_ents_to_finite_element_by_dim
virtual MoFEMErrorCode add_ents_to_finite_element_by_dim(const EntityHandle entities, const int dim, const std::string &name, const bool recursive=true)=0
add entities to finite element
MoFEM::FatPrismElementForcesAndSourcesCore
FatPrism finite element.
Definition: FatPrismElementForcesAndSourcesCore.hpp:31
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator
Definition: VolumeElementForcesAndSourcesCore.hpp:108
FTensor::Index< 'i', 3 >
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1536
convert.n
n
Definition: convert.py:82
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1518
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
MF_ZERO
@ MF_ZERO
Definition: definitions.h:98
MoFEM::EntitiesFieldData::EntData::getFTensor1DiffN
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
Definition: EntitiesFieldData.cpp:526
FTensor::Tensor0
Definition: Tensor0.hpp:16
MoFEM::EntitiesFieldData::EntData::getN
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
Definition: EntitiesFieldData.hpp:1305
OpCalculateStrain
Definition: HookeElement.hpp:119
MoFEM::DMoFEMGetInterfacePtr
PetscErrorCode DMoFEMGetInterfacePtr(DM dm, MoFEM::Interface **m_field_ptr)
Get pointer to MoFEM::Interface.
Definition: DMMoFEM.cpp:414
MAT_TO_DDG
#define MAT_TO_DDG(SM)
Definition: HookeElement.hpp:142
BLOCKSET
@ BLOCKSET
Definition: definitions.h:148
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
FTensor::Ddg
Definition: Ddg_value.hpp:7
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
OpAleLhs_dx_dX
Definition: HookeElement.hpp:449
MoFEM::createVectorMPI
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
Definition: PetscSmartObj.hpp:202
MoFEM::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
addElasticElement
static MoFEMErrorCode addElasticElement(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData >> &block_sets_ptr, const std::string element_name, const std::string x_field, const std::string X_field, const bool ale)
Definition: HookeElement.cpp:533
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
OpCalculateStress
Definition: HookeElement.hpp:153
setOperators
static MoFEMErrorCode setOperators(boost::shared_ptr< ForcesAndSourcesCore > fe_lhs_ptr, boost::shared_ptr< ForcesAndSourcesCore > fe_rhs_ptr, boost::shared_ptr< map< int, BlockData >> block_sets_ptr, const std::string x_field, const std::string X_field, const bool ale, const bool field_disp, const EntityType type=MBTET, boost::shared_ptr< DataAtIntegrationPts > data_at_pts=nullptr)
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
OpAssemble
Definition: HookeElement.hpp:340
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
OpAleLhs_dx_dx
Definition: HookeElement.hpp:434
MoFEM::SmartPetscObj< Vec >
OpAleLhs_dX_dX
Definition: HookeElement.hpp:521
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:586
OpRhs_dx
Definition: HookeElement.hpp:401
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
F
@ F
Definition: free_surface.cpp:394
OpAleLhs_dX_dx
Definition: HookeElement.hpp:547