Operators and data structures for linear elastic analysis.
In other words spatial deformation is small but topological changes large.
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_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 determinantTensor3by3(t_H, t_detH);
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,
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];
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) {
if (type == MBPRISM) {
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 {
if (type == MBPRISM) {
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) {
if (type == MBPRISM) {
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 {
if (type == MBPRISM) {
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) {
if (type == MBPRISM) {
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 {
if (type == MBPRISM) {
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 rho = getFTensor0FromVec(*rhoAtGaussPtsPtr);
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 rho = getFTensor0FromVec(*rhoAtGaussPtsPtr);
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;
auto rho = getFTensor0FromVec(*rhoAtGaussPtsPtr);
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;
}
}
}
#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 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_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
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_row(const std::string &fe_name, const std::string name_row)=0
set field row 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 c
speed of light (cm/ns)
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 createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
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.
intrusive_ptr for managing petsc objects