Operators and data structures for linear elastic analysis.
Operators and data structures for linear elastic analysisSee as well header file HookeElement.hpp
In other words spatial deformation is small but topological changes large.
#include <MoFEM.hpp>
HookeElement::OpCalculateStrainAle::OpCalculateStrainAle(
const std::string row_field, const std::string col_field,
boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
dataAtPts(data_at_pts) {
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
MoFEMErrorCode HookeElement::OpCalculateStrainAle::doWork(
int row_side,
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_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) {
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,
dataAtPts(data_at_pts), ghostVec(ghost_vec, true) {
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
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()) {
double vol = getVolume();
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) {
if (det_H.size()) {
}
++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)
dataAtPts(data_at_pts) {
std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
}
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)
dataAtPts(data_at_pts) {}
MoFEMErrorCode HookeElement::OpAssemble::doWork(
int row_side,
int col_side,
if (!nbRows)
if (!nbCols)
K.resize(nbRows, nbCols,
false);
nbIntegrationPts = getGaussPts().size2();
if (row_side == col_side && row_type == col_type) {
isDiag = true;
} else {
isDiag = false;
}
CHKERR iNtegrate(row_data, col_data);
CHKERR aSsemble(row_data, col_data);
}
if (!nbRows)
nF.resize(nbRows, false);
nF.clear();
nbIntegrationPts = getGaussPts().size2();
}
};
};
const int *row_indices = &*row_data.
getIndices().data().begin();
const int *col_indices = &*col_data.
getIndices().data().begin();
auto &data = *dataAtPts;
if (!data.forcesOnlyOnEntitiesRow.empty()) {
rowIndices.resize(nbRows, false);
row_indices = &rowIndices[0];
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);
col_indices = &colIndices[0];
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;
&*
K.data().begin(), ADD_VALUES);
if (!isDiag && sYmm) {
transK.resize(
K.size2(),
K.size1(),
false);
noalias(transK) = trans(
K);
&*transK.data().begin(), ADD_VALUES);
}
}
const int *row_indices = &*row_data.
getIndices().data().begin();
auto &data = *dataAtPts;
if (!data.forcesOnlyOnEntitiesRow.empty()) {
rowIndices.resize(nbRows, false);
row_indices = &rowIndices[0];
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;
}
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) {}
&
v(
r + 0), &
v(
r + 1), &
v(
r + 2));
};
double vol = getVolume();
auto t_w = getFTensor0IntegrationWeight();
const int row_nb_base_fun = row_data.
getN().size2();
auto t_cauchy_stress =
getFTensor2SymmetricFromMat<3>(*dataAtPts->cauchyStressMat);
for (int gg = 0; gg != nbIntegrationPts; ++gg) {
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) {}
&
v(
r + 0), &
v(
r + 1), &
v(
r + 2));
};
double vol = getVolume();
auto t_w = getFTensor0IntegrationWeight();
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);
for (int gg = 0; gg != nbIntegrationPts; ++gg) {
double a = t_w * vol * det_H[gg];
auto t_nf = get_tensor1(nF, 0);
int rr = 0;
for (; rr != nbRows / 3; ++rr) {
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) {}
&
v(
r + 0), &
v(
r + 1), &
v(
r + 2));
};
double vol = getVolume();
auto t_w = getFTensor0IntegrationWeight();
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);
for (int gg = 0; gg != nbIntegrationPts; ++gg) {
double a = t_w * vol * det_H[gg];
auto t_nf = get_tensor1(nF, 0);
int rr = 0;
for (; rr != nbRows / 3; ++rr) {
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;
}
}
boost::shared_ptr<map<int, BlockData>> &block_sets_ptr) {
if (!block_sets_ptr)
"Pointer to block of sets is null");
CHKERR it->getAttributeDataStructure(mydata);
int id = it->getMeshsetId();
auto &block_data = (*block_sets_ptr)[id];
EntityHandle meshset = it->getMeshset();
block_data.tEts, true);
block_data.iD = id;
block_data.E = mydata.
data.Young;
block_data.PoissonRatio = mydata.
data.Poisson;
}
}
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)
"Pointer to block of sets is null");
if (ale) {
}
}
for (
auto &
m : (*block_sets_ptr)) {
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,
boost::shared_ptr<DataAtIntegrationPts> data_at_pts) {
if (!block_sets_ptr)
"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) {
boost::shared_ptr<MatrixDouble> inv_jac_ptr(
new MatrixDouble);
fe_lhs_ptr->getOpPtrVector().push_back(
fe_lhs_ptr->getOpPtrVector().push_back(
}
fe_lhs_ptr->getOpPtrVector().push_back(
x_field, x_field, block_sets_ptr, data_at_pts));
fe_lhs_ptr->getOpPtrVector().push_back(
} else {
boost::shared_ptr<MatrixDouble> inv_jac_ptr(
new MatrixDouble);
fe_lhs_ptr->getOpPtrVector().push_back(
fe_lhs_ptr->getOpPtrVector().push_back(
}
fe_lhs_ptr->getOpPtrVector().push_back(
fe_lhs_ptr->getOpPtrVector().push_back(
x_field, x_field, block_sets_ptr, data_at_pts));
fe_lhs_ptr->getOpPtrVector().push_back(
fe_lhs_ptr->getOpPtrVector().push_back(
fe_lhs_ptr->getOpPtrVector().push_back(
fe_lhs_ptr->getOpPtrVector().push_back(
fe_lhs_ptr->getOpPtrVector().push_back(
fe_lhs_ptr->getOpPtrVector().push_back(
fe_lhs_ptr->getOpPtrVector().push_back(
fe_lhs_ptr->getOpPtrVector().push_back(
fe_lhs_ptr->getOpPtrVector().push_back(
fe_lhs_ptr->getOpPtrVector().push_back(
}
}
if (fe_rhs_ptr) {
if (ale == PETSC_FALSE) {
boost::shared_ptr<MatrixDouble> inv_jac_ptr(
new MatrixDouble);
fe_rhs_ptr->getOpPtrVector().push_back(
fe_rhs_ptr->getOpPtrVector().push_back(
}
fe_rhs_ptr->getOpPtrVector().push_back(
fe_rhs_ptr->getOpPtrVector().push_back(
block_sets_ptr, data_at_pts));
if (field_disp) {
fe_rhs_ptr->getOpPtrVector().push_back(
} else {
fe_rhs_ptr->getOpPtrVector().push_back(
}
fe_rhs_ptr->getOpPtrVector().push_back(
fe_rhs_ptr->getOpPtrVector().push_back(
new OpRhs_dx(x_field, x_field, data_at_pts));
} else {
boost::shared_ptr<MatrixDouble> inv_jac_ptr(
new MatrixDouble);
fe_rhs_ptr->getOpPtrVector().push_back(
fe_rhs_ptr->getOpPtrVector().push_back(
}
fe_rhs_ptr->getOpPtrVector().push_back(
fe_rhs_ptr->getOpPtrVector().push_back(
block_sets_ptr, data_at_pts));
fe_rhs_ptr->getOpPtrVector().push_back(
fe_rhs_ptr->getOpPtrVector().push_back(
fe_rhs_ptr->getOpPtrVector().push_back(
fe_rhs_ptr->getOpPtrVector().push_back(
fe_rhs_ptr->getOpPtrVector().push_back(
fe_rhs_ptr->getOpPtrVector().push_back(
fe_rhs_ptr->getOpPtrVector().push_back(
}
}
}
DM dm, boost::shared_ptr<map<int, BlockData>> block_sets_ptr,
const std::string x_field, const std::string X_field, const bool ale,
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) {
};
"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(
auto push_ops = [&](boost::shared_ptr<ForcesAndSourcesCore> fe_ptr,
boost::shared_ptr<MatrixDouble> inv_jac_ptr(
new MatrixDouble);
if (ale == PETSC_FALSE) {
fe_ptr->getOpPtrVector().push_back(
fe_ptr->getOpPtrVector().push_back(
}
fe_ptr->getOpPtrVector().push_back(
x_field, x_field, block_sets_ptr, data_at_pts));
if (field_disp) {
fe_ptr->getOpPtrVector().push_back(
} else {
fe_ptr->getOpPtrVector().push_back(
}
fe_ptr->getOpPtrVector().push_back(
fe_ptr->getOpPtrVector().push_back(
} else {
fe_ptr->getOpPtrVector().push_back(
fe_ptr->getOpPtrVector().push_back(
}
fe_ptr->getOpPtrVector().push_back(
x_field, x_field, block_sets_ptr, data_at_pts));
fe_ptr->getOpPtrVector().push_back(
fe_ptr->getOpPtrVector().push_back(
fe_ptr->getOpPtrVector().push_back(
fe_ptr->getOpPtrVector().push_back(
}
};
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 VecAssemblyBegin(v_energy);
CHKERR VecAssemblyEnd(v_energy);
}
&
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),
};
double vol = getVolume();
auto t_w = getFTensor0IntegrationWeight();
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();
for (int gg = 0; gg != nbIntegrationPts; ++gg) {
double a = t_w * vol * det_H[gg];
int rr = 0;
for (; rr != nbRows / 3; ++rr) {
auto t_m = get_tensor2(
K, 3 * rr, 0);
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);
for (int cc = 0; cc != nbCols / 3; ++cc) {
t_m(
i,
k) += t_row_stress_dx(
i,
k,
l) * t_col_diff_base(
l);
++t_col_diff_base;
++t_m;
}
++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,
&
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),
};
double vol = getVolume();
auto t_w = getFTensor0IntegrationWeight();
const int row_nb_base_fun = row_data.
getN().size2();
auto t_grad_rho = getFTensor1FromMat<3>(*rhoGradAtGaussPtsPtr);
auto t_eshelby_stress =
getFTensor2FromMat<3, 3>(*dataAtPts->eshelbyStressMat);
auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
auto &det_H = *dataAtPts->detHVec;
for (int gg = 0; gg != nbIntegrationPts; ++gg) {
double a = t_w * vol * det_H[gg];
const double stress_dho_coef = (rhoN /
rho);
int rr = 0;
for (; rr != nbRows / 3; ++rr) {
auto t_m = get_tensor2(
K, 3 * rr, 0);
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);
for (int cc = 0; cc != nbCols / 3; ++cc) {
t_row_stress(
i) * stress_dho_coef * t_grad_rho(
k) * t_col_base;
++t_col_base;
++t_m;
}
++t_row_diff_base;
}
for (; rr != row_nb_base_fun; ++rr)
++t_row_diff_base;
++t_w;
++t_eshelby_stress;
++t_invH;
++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,
&
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),
};
double vol = getVolume();
auto t_w = getFTensor0IntegrationWeight();
const int row_nb_base_fun = row_data.
getN().size2();
auto t_grad_rho = getFTensor1FromMat<3>(*rhoGradAtGaussPtsPtr);
auto t_cauchy_stress =
getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
auto &det_H = *dataAtPts->detHVec;
for (int gg = 0; gg != nbIntegrationPts; ++gg) {
double a = t_w * vol * det_H[gg];
const double stress_dho_coef = (rhoN /
rho);
int rr = 0;
for (; rr != nbRows / 3; ++rr) {
auto t_m = get_tensor2(
K, 3 * rr, 0);
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);
for (int cc = 0; cc != nbCols / 3; ++cc) {
t_row_stress(
i) * stress_dho_coef * t_grad_rho(
k) * t_col_base;
++t_col_base;
++t_m;
}
++t_row_diff_base;
}
for (; rr != row_nb_base_fun; ++rr)
++t_row_diff_base;
++t_w;
++t_cauchy_stress;
++t_invH;
++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)
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(
if (!rhoAtGaussPtsPtr)
SETERRQ(PETSC_COMM_SELF, 1, "Calculate density with MWLS first.");
for (
auto &
m : (*blockSetsPtr)) {
if (
m.second.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
continue;
}
const int nb_integration_pts = getGaussPts().size2();
dataAtPts->stiffnessMat->resize(36, nb_integration_pts, false);
const double young =
m.second.E;
const double poisson =
m.second.PoissonRatio;
const double coefficient = young / ((1 + poisson) * (1 - 2 * poisson));
for (int gg = 0; gg != nb_integration_pts; ++gg) {
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;
t_D(
i,
j,
k,
l) *= coefficient * pow(
rho / rHo0, rhoN);
++t_D;
}
}
}
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)
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)
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
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.
PetscErrorCode DMoFEMGetInterfacePtr(DM dm, MoFEM::Interface **m_field_ptr)
Get pointer to MoFEM::Interface.
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.
UBlasMatrix< double > MatrixDouble
UBlasVector< double > VectorDouble
implementation of Data Operators for Forces and Sources
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.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
constexpr auto VecSetValues
constexpr auto MatSetValues
const double r
rate factor
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.