v0.13.1
HookeElement.cpp

Operators and data structures for linear elastic analysis.

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 = createSmartVectorMPI(m_field_ptr->get_comm(), PETSC_DECIDE, 1);
boost::shared_ptr<DataAtIntegrationPts> data_at_pts(
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);
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,
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;
}
}
}
static Index< 'o', 3 > o
static Index< 'K', 3 > K
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)
#define MAT_TO_DDG(SM)
static MoFEMErrorCode setBlocks(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData > > &block_sets_ptr)
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)
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)
constexpr double a
@ MF_ZERO
Definition: definitions.h:98
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
Definition: definitions.h:159
@ BLOCKSET
Definition: definitions.h:148
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:533
PetscErrorCode DMoFEMGetInterfacePtr(DM dm, MoFEM::Interface **m_field_ptr)
Get pointer to MoFEM::Interface.
Definition: DMMMoFEM.cpp:361
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
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
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
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
virtual bool check_field(const std::string &name) const =0
check if field is in database
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
auto createSmartVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
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:1204
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1193
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
const double r
rate factor
double rho
Definition: plastic.cpp:89
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Elastic material data structure.
Calculate inverse of jacobian for face element.
Transform local reference derivatives of shape functions to global derivatives.