#include <boost/math/constants/constants.hpp>
#ifdef PYTHON_SDF
#include <boost/python.hpp>
#include <boost/python/def.hpp>
#include <boost/python/numpy.hpp>
namespace bp = boost::python;
namespace np = boost::python::numpy;
#pragma message "With PYTHON_SDF"
#else
#pragma message "Without PYTHON_SDF"
#endif
struct VolUserDataOperatorStabAssembly
};
struct FaceUserDataOperatorStabAssembly
};
}
const std::string block_name) {
std::regex((boost::format("%s(.*)") % block_name).str())
);
for (auto bc : bcs) {
bc->getMeshsetIdEntitiesByDimension(m_field.
get_moab(), 2, faces,
true),
"get meshset ents");
}
};
template <>
return MatSetValues<AssemblyTypeSelector<SCHUR>>(
op_ptr->
getKSPB(), row_data, col_data,
m, ADD_VALUES);
};
PetscBool EshelbianCore::noStretch = PETSC_FALSE;
double EshelbianCore::exponentBase = exp(1);
EshelbianCore::d_f_log;
EshelbianCore::dd_f_log;
boost::function<
double(
const double)> EshelbianCore::inv_f =
EshelbianCore::inv_f_log;
boost::function<
double(
const double)> EshelbianCore::inv_d_f =
EshelbianCore::inv_d_f_log;
boost::function<
double(
const double)> EshelbianCore::inv_dd_f =
EshelbianCore::inv_dd_f_log;
EshelbianCore::query_interface(boost::typeindex::type_index type_index,
return 0;
}
if (evalRhs)
if (evalLhs)
}
: mField(m_field), piolaStress("P"), eshelbyStress("S"),
spatialL2Disp("wL2"), spatialH1Disp("wH1"), contactDisp("contactDisp"),
materialL2Disp("W"), stretchTensor("u"), rotAxis("omega"),
materialGradient("G"), tauField("TAU"), lambdaField("LAMBDA"),
bubbleField("BUBBLE"), elementVolumeName("EP"),
naturalBcElement("NATURAL_BC"), essentialBcElement("ESSENTIAL_BC"),
skinElement("SKIN_ELEMENT"), contactElement("CONTACT_ELEMENT") {
CHKERRABORT(PETSC_COMM_WORLD,
ierr);
}
EshelbianCore::~EshelbianCore() = default;
const char *list_rots[] = {"small", "moderate", "large", "no_h1"};
PetscInt choice_rot = EshelbianCore::rotSelector;
PetscInt choice_grad = EshelbianCore::gradApperoximator;
const char *list_stretches[] = {"linear", "log"};
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Eshelbian plasticity",
"none");
CHKERR PetscOptionsInt(
"-space_order",
"approximation oder for space",
"",
spaceOrder, &spaceOrder, PETSC_NULL);
CHKERR PetscOptionsInt(
"-space_h1_order",
"approximation oder for space",
"",
spaceH1Order, &spaceH1Order, PETSC_NULL);
CHKERR PetscOptionsInt(
"-material_order",
"approximation oder for material",
"", materialOrder, &materialOrder, PETSC_NULL);
CHKERR PetscOptionsScalar(
"-viscosity_alpha_u",
"viscosity",
"", alphaU,
&alphaU, PETSC_NULL);
CHKERR PetscOptionsScalar(
"-viscosity_alpha_w",
"viscosity",
"", alphaW,
&alphaW, PETSC_NULL);
CHKERR PetscOptionsScalar(
"-density_alpha_rho",
"density",
"", alphaRho,
&alphaRho, PETSC_NULL);
CHKERR PetscOptionsScalar(
"-preconditioner_eps",
"preconditioner_eps",
"",
precEps, &precEps, PETSC_NULL);
CHKERR PetscOptionsScalar(
"-preconditioner_eps_omega",
"preconditioner_eps",
"", precEpsOmega, &precEpsOmega, PETSC_NULL);
CHKERR PetscOptionsScalar(
"-preconditioner_eps_w",
"preconditioner_eps",
"",
precEpsW, &precEpsW, PETSC_NULL);
CHKERR PetscOptionsScalar(
"-preconditioner_eps_contact_disp",
"preconditioner_eps", "", precEpsContactDisp,
&precEpsContactDisp, PETSC_NULL);
CHKERR PetscOptionsEList(
"-rotations",
"rotations",
"", list_rots,
LARGE_ROT + 1, list_rots[choice_rot], &choice_rot,
PETSC_NULL);
CHKERR PetscOptionsEList(
"-grad",
"gradient of defamation approximate",
"",
list_rots[choice_grad], &choice_grad, PETSC_NULL);
CHKERR PetscOptionsScalar(
"-exponent_base",
"exponent_base",
"", exponentBase,
&EshelbianCore::exponentBase, PETSC_NULL);
list_stretches[choice_stretch], &choice_stretch, PETSC_NULL);
CHKERR PetscOptionsBool(
"-no_stretch",
"do not solve for stretch",
"",
noStretch, &noStretch, PETSC_NULL);
ierr = PetscOptionsEnd();
EshelbianCore::rotSelector =
static_cast<RotSelector>(choice_rot);
EshelbianCore::gradApperoximator =
static_cast<RotSelector>(choice_grad);
EshelbianCore::stretchSelector =
static_cast<StretchSelector>(choice_stretch);
switch (EshelbianCore::stretchSelector) {
EshelbianCore::inv_f = inv_f_linear;
EshelbianCore::inv_d_f = inv_d_f_linear;
EshelbianCore::inv_dd_f = inv_dd_f_linear;
break;
EshelbianCore::inv_f = inv_f_log;
EshelbianCore::inv_d_f = inv_d_f_log;
EshelbianCore::inv_dd_f = inv_dd_f_log;
break;
default:
break;
};
precEpsOmega += precEps;
precEpsW += precEps;
MOFEM_LOG(
"EP", Sev::inform) <<
"spaceOrder " << spaceOrder;
MOFEM_LOG(
"EP", Sev::inform) <<
"spaceH1Order " << spaceH1Order;
MOFEM_LOG(
"EP", Sev::inform) <<
"materialOrder " << materialOrder;
MOFEM_LOG(
"EP", Sev::inform) <<
"alphaU " << alphaU;
MOFEM_LOG(
"EP", Sev::inform) <<
"alphaW " << alphaW;
MOFEM_LOG(
"EP", Sev::inform) <<
"alphaRho " << alphaRho;
MOFEM_LOG(
"EP", Sev::inform) <<
"precEps " << precEps;
MOFEM_LOG(
"EP", Sev::inform) <<
"precEpsOmega " << precEpsOmega;
MOFEM_LOG(
"EP", Sev::inform) <<
"precEpsW " << precEpsW;
<< "Rotations " << list_rots[EshelbianCore::rotSelector];
MOFEM_LOG(
"EP", Sev::inform) <<
"Gradient of deformation "
<< list_rots[EshelbianCore::gradApperoximator];
if (exponentBase != exp(1))
<< "Base exponent " << EshelbianCore::exponentBase;
else
MOFEM_LOG(
"EP", Sev::inform) <<
"Base exponent e";
MOFEM_LOG(
"EP", Sev::inform) <<
"Stretch " << list_stretches[choice_stretch];
if (spaceH1Order == -1)
spaceH1Order = spaceOrder;
}
CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets);
Skinner skin(&mField.get_moab());
CHKERR skin.find_skin(0, tets,
false, tets_skin_part);
ParallelComm *pcomm =
CHKERR pcomm->filter_pstatus(tets_skin_part,
PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &tets_skin);
if (bcSpatialDispVecPtr)
for (
auto &
v : *bcSpatialDispVecPtr) {
tets_skin = subtract(tets_skin,
v.faces);
}
if (bcSpatialRotationVecPtr)
for (
auto &
v : *bcSpatialRotationVecPtr) {
tets_skin = subtract(tets_skin,
v.faces);
}
if (bcSpatialTraction)
for (
auto &
v : *bcSpatialTraction) {
tets_skin = subtract(tets_skin,
v.faces);
}
auto subtract_faces_where_displacements_are_applied =
[&](const std::string block_name) {
tets_skin = subtract(tets_skin, contact_range);
};
CHKERR subtract_faces_where_displacements_are_applied(
"CONTACT");
moab::Interface::UNION);
Range faces_not_on_the_skin = subtract(faces, tets_skin);
const int dim) {
auto field_ptr = mField.get_field_structure(
field_name);
auto meshset = field_ptr->getMeshset();
CHKERR mField.get_moab().add_entities(meshset, tets);
CHKERR mField.get_moab().add_entities(meshset, faces_not_on_the_skin);
};
auto add_l2_field = [
this, meshset](
const std::string
field_name,
const int order,
const int dim) {
};
auto add_h1_field = [
this, meshset](
const std::string
field_name,
const int order,
const int dim) {
};
auto add_l2_field_by_range = [
this](
const std::string
field_name,
const int order,
const int dim,
const int field_dim,
Range &&
r) {
};
auto add_bubble_field = [
this, meshset](
const std::string
field_name,
const int order,
const int dim) {
auto field_ptr = mField.get_field_structure(
field_name);
auto field_order_table =
const_cast<Field *
>(field_ptr)->getFieldOrderTable();
auto get_cgg_bubble_order_zero = [](int p) { return 0; };
auto get_cgg_bubble_order_tet = [](int p) {
};
field_order_table[MBVERTEX] = get_cgg_bubble_order_zero;
field_order_table[MBEDGE] = get_cgg_bubble_order_zero;
field_order_table[MBTRI] = get_cgg_bubble_order_zero;
field_order_table[MBTET] = get_cgg_bubble_order_tet;
};
CHKERR add_hdiv_field(piolaStress, spaceOrder, 3);
CHKERR add_bubble_field(bubbleField, spaceOrder, 1);
CHKERR add_l2_field(spatialL2Disp, spaceOrder - 1, 3);
CHKERR add_l2_field(rotAxis, spaceOrder - 1, 3);
CHKERR add_l2_field(stretchTensor, spaceOrder, 6);
CHKERR add_l2_field_by_range(contactDisp, spaceOrder - 1, 2, 3,
CHKERR add_h1_field(spatialH1Disp, spaceH1Order, 3);
}
EshelbianCore::addVolumeFiniteElement(
const EntityHandle meshset) {
auto add_field_to_fe = [this](const std::string fe,
};
if (!mField.check_finite_element(elementVolumeName)) {
CHKERR mField.add_ents_to_finite_element_by_type(meshset, MBTET,
elementVolumeName);
CHKERR add_field_to_fe(elementVolumeName, piolaStress);
CHKERR add_field_to_fe(elementVolumeName, bubbleField);
CHKERR add_field_to_fe(elementVolumeName, stretchTensor);
CHKERR add_field_to_fe(elementVolumeName, rotAxis);
CHKERR add_field_to_fe(elementVolumeName, spatialL2Disp);
CHKERR add_field_to_fe(elementVolumeName, spatialH1Disp);
CHKERR add_field_to_fe(elementVolumeName, contactDisp);
}
CHKERR mField.build_finite_elements(elementVolumeName);
}
EshelbianCore::addBoundaryFiniteElement(
const EntityHandle meshset) {
auto add_field_to_fe = [this](const std::string fe,
};
Range natural_bc_elements;
if (bcSpatialDispVecPtr) {
for (
auto &
v : *bcSpatialDispVecPtr) {
natural_bc_elements.merge(
v.faces);
}
}
if (bcSpatialRotationVecPtr) {
for (
auto &
v : *bcSpatialRotationVecPtr) {
natural_bc_elements.merge(
v.faces);
}
}
Range essentail_bc_elements;
if (bcSpatialTraction) {
for (
auto &
v : *bcSpatialTraction) {
essentail_bc_elements.merge(
v.faces);
}
}
CHKERR mField.add_ents_to_finite_element_by_type(natural_bc_elements, MBTRI,
naturalBcElement);
CHKERR add_field_to_fe(naturalBcElement, piolaStress);
CHKERR mField.build_finite_elements(naturalBcElement);
CHKERR mField.add_ents_to_finite_element_by_type(essentail_bc_elements, MBTRI,
essentialBcElement);
CHKERR add_field_to_fe(essentialBcElement, piolaStress);
CHKERR mField.build_finite_elements(essentialBcElement);
auto get_skin = [&]() {
CHKERR mField.get_moab().get_entities_by_dimension(0,
SPACE_DIM, body_ents);
Skinner skin(&mField.get_moab());
CHKERR skin.find_skin(0, body_ents,
false, skin_ents);
return skin_ents;
};
auto filter_true_skin = [&](auto &&skin) {
ParallelComm *pcomm =
CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &boundary_ents);
return boundary_ents;
};
auto skin = filter_true_skin(get_skin());
CHKERR mField.add_ents_to_finite_element_by_type(skin, MBTRI, skinElement);
CHKERR add_field_to_fe(skinElement, piolaStress);
CHKERR add_field_to_fe(skinElement, spatialH1Disp);
CHKERR add_field_to_fe(skinElement, contactDisp);
CHKERR mField.build_finite_elements(skinElement);
MOFEM_LOG(
"EP", Sev::inform) <<
"Contact elements " << contact_range.size();
CHKERR mField.add_ents_to_finite_element_by_type(contact_range, MBTRI,
contactElement);
CHKERR add_field_to_fe(contactElement, piolaStress);
CHKERR add_field_to_fe(contactElement, spatialH1Disp);
CHKERR add_field_to_fe(contactElement, contactDisp);
CHKERR mField.build_finite_elements(contactElement);
}
dM =
createDM(mField.get_comm(),
"DMMOFEM");
mField.getInterface<
ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
mField.getInterface<
ProblemsManager>()->buildProblemFromFields = PETSC_FALSE;
auto remove_dofs_on_essential_spatial_stress_boundary =
[&](const std::string prb_name) {
all_three.merge((*bcSpatialFreeTraction)[
d]);
all_three = intersect(all_three, (*bcSpatialFreeTraction)[
d]);
prb_name, piolaStress, all_three, 0,
SPACE_DIM, 0,
};
CHKERR remove_dofs_on_essential_spatial_stress_boundary(
"ESHELBY_PLASTICITY");
dmElastic =
createDM(mField.get_comm(),
"DMMOFEM");
if (!noStretch) {
}
auto set_zero_block = [&]() {
"ELASTIC_PROBLEM", spatialL2Disp, stretchTensor);
"ELASTIC_PROBLEM", stretchTensor, spatialL2Disp);
"ELASTIC_PROBLEM", spatialL2Disp, rotAxis);
"ELASTIC_PROBLEM", rotAxis, spatialL2Disp);
"ELASTIC_PROBLEM", spatialL2Disp, bubbleField);
"ELASTIC_PROBLEM", bubbleField, spatialL2Disp);
"ELASTIC_PROBLEM", bubbleField, bubbleField);
"ELASTIC_PROBLEM", bubbleField, piolaStress);
"ELASTIC_PROBLEM", piolaStress, bubbleField);
};
auto set_section = [&]() {
PetscSection section;
§ion);
CHKERR DMSetSection(dmElastic, section);
CHKERR DMSetGlobalSection(dmElastic, section);
CHKERR PetscSectionDestroy(§ion);
};
dmPrjSpatial =
createDM(mField.get_comm(),
"DMMOFEM");
->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
"PROJECT_SPATIAL", spatialH1Disp, true, false);
}
BcDisp::BcDisp(std::string name, std::vector<double> &attr,
Range &faces)
: blockName(name), faces(faces) {
vals.resize(3, false);
flags.resize(3, false);
for (int ii = 0; ii != 3; ++ii) {
vals[ii] = attr[ii];
flags[ii] = static_cast<int>(attr[ii + 3]);
}
MOFEM_LOG(
"EP", Sev::inform) <<
"Add BCDisp " << name;
<< "Add BCDisp vals " << vals[0] << " " << vals[1] << " " << vals[2];
<< "Add BCDisp flags " << flags[0] << " " << flags[1] << " " << flags[2];
MOFEM_LOG(
"EP", Sev::inform) <<
"Add BCDisp nb. of faces " << faces.size();
}
BcRot::BcRot(std::string name, std::vector<double> &attr,
Range &faces)
: blockName(name), faces(faces) {
vals.resize(3, false);
for (int ii = 0; ii != 3; ++ii) {
vals[ii] = attr[ii];
}
theta = attr[3];
}
TractionBc::TractionBc(std::string name, std::vector<double> &attr,
: blockName(name), faces(faces) {
vals.resize(3, false);
flags.resize(3, false);
for (int ii = 0; ii != 3; ++ii) {
vals[ii] = attr[ii];
flags[ii] = static_cast<int>(attr[ii + 3]);
}
MOFEM_LOG(
"EP", Sev::inform) <<
"Add BCForce " << name;
<< "Add BCForce vals " << vals[0] << " " << vals[1] << " " << vals[2];
<< "Add BCForce flags " << flags[0] << " " << flags[1] << " " << flags[2];
MOFEM_LOG(
"EP", Sev::inform) <<
"Add BCForce nb. of faces " << faces.size();
}
EshelbianCore::getTractionFreeBc(
const EntityHandle meshset,
boost::shared_ptr<TractionFreeBc> &bc_ptr,
const std::string contact_set_name) {
CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets);
Skinner skin(&mField.get_moab());
CHKERR skin.find_skin(0, tets,
false, tets_skin_part);
ParallelComm *pcomm =
CHKERR pcomm->filter_pstatus(tets_skin_part,
PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &tets_skin);
bc_ptr->resize(3);
for (
int dd = 0;
dd != 3; ++
dd)
(*bc_ptr)[
dd] = tets_skin;
if (bcSpatialDispVecPtr)
for (
auto &
v : *bcSpatialDispVecPtr) {
(*bc_ptr)[0] = subtract((*bc_ptr)[0],
v.faces);
(*bc_ptr)[1] = subtract((*bc_ptr)[1],
v.faces);
(*bc_ptr)[2] = subtract((*bc_ptr)[2],
v.faces);
}
if (bcSpatialRotationVecPtr)
for (
auto &
v : *bcSpatialRotationVecPtr) {
(*bc_ptr)[0] = subtract((*bc_ptr)[0],
v.faces);
(*bc_ptr)[1] = subtract((*bc_ptr)[1],
v.faces);
(*bc_ptr)[2] = subtract((*bc_ptr)[2],
v.faces);
}
std::regex((boost::format("%s(.*)") % contact_set_name).str()))) {
CHKERR m->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
true);
(*bc_ptr)[0] = subtract((*bc_ptr)[0], faces);
(*bc_ptr)[1] = subtract((*bc_ptr)[1], faces);
(*bc_ptr)[2] = subtract((*bc_ptr)[2], faces);
}
}
int operator()(
int p_row,
int p_col,
int p_data)
const {
return p_data + p_data + (p_data - 1);
else
return 2 * p_data;
}
};
int operator()(
int p_row,
int p_col,
int p_data)
const {
return 2 * p_data; }
};
CGGUserPolynomialBase() {}
~CGGUserPolynomialBase() {}
*iface = const_cast<CGGUserPolynomialBase *>(this);
return 0;
}
boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
int nb_gauss_pts = pts.size2();
if (!nb_gauss_pts) {
}
if (pts.size1() < 3) {
"Wrong dimension of pts, should be at least 3 rows with "
"coordinates");
}
switch (cTx->sPace) {
CHKERR getValueHdivForCGGBubble(pts);
break;
default:
}
}
private:
boost::tuple<int, int, MatrixDouble> cachePhi = {0, 0,
MatrixDouble()};
"Wrong base, should be USER_BASE");
}
const int nb_gauss_pts = pts.size2();
if (cachePhi.get<0>() ==
order && cachePhi.get<1>() == nb_gauss_pts) {
auto &mat = cachePhi.get<2>();
phi.resize(mat.size1(), mat.size2(),
false);
} else {
shapeFun.resize(nb_gauss_pts, 4, false);
&pts(2, 0), nb_gauss_pts);
double diff_shape_fun[12];
phi.resize(nb_gauss_pts, 9 * nb_base_functions,
false);
nb_gauss_pts);
cachePhi.get<0>() =
order;
cachePhi.get<1>() = nb_gauss_pts;
cachePhi.get<2>().resize(
phi.size1(),
phi.size2(),
false);
noalias(cachePhi.get<2>()) =
phi;
}
}
};
const int tag, const bool do_rhs, const bool do_lhs,
boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe) {
fe = boost::make_shared<VolumeElementForcesAndSourcesCore>(mField);
fe->getUserPolynomialBase() =
boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
{HDIV, H1, L2});
if (!dataAtPts) {
dataAtPts =
boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
dataAtPts->physicsPtr = physicalEquations;
}
piolaStress, dataAtPts->getApproxPAtPts()));
bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
piolaStress, dataAtPts->getDivPAtPts()));
if (noStretch) {
fe->getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
tag, do_rhs, do_lhs, dataAtPts, physicalEquations));
fe->getOpPtrVector().push_back(new OpCalculateStretchFromStress(dataAtPts));
} else {
fe->getOpPtrVector().push_back(
stretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
}
rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
spatialL2Disp, dataAtPts->getSmallWL2AtPts(), MBTET));
spatialL2Disp, dataAtPts->getSmallWL2DotAtPts(), MBTET));
if (!noStretch) {
fe->getOpPtrVector().push_back(
stretchTensor, dataAtPts->getLogStretchDotTensorAtPts(), MBTET));
}
rotAxis, dataAtPts->getRotAxisDotAtPts(), MBTET));
if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
spatialL2Disp, dataAtPts->getSmallWL2DotDotAtPts(), MBTET));
}
spatialH1Disp, dataAtPts->getSmallWH1AtPts()));
spatialH1Disp, dataAtPts->getSmallWGradH1AtPts()));
fe->getOpPtrVector().push_back(
new OpCalculateRotationAndSpatialGradient(dataAtPts));
if (!noStretch) {
fe->getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
tag, do_rhs, do_lhs, dataAtPts, physicalEquations));
}
}
const int tag, const bool add_elastic, const bool add_material,
boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe_rhs,
boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe_lhs) {
auto time_scale = boost::make_shared<TimeScale>();
CHKERR setBaseVolumeElementOps(tag,
true,
false, fe_rhs);
if (add_elastic) {
fe_rhs->getOpPtrVector().push_back(
new OpSpatialEquilibrium(spatialL2Disp, dataAtPts, alphaW, alphaRho));
fe_rhs->getOpPtrVector().push_back(
new OpSpatialRotation(rotAxis, dataAtPts));
if (!noStretch) {
fe_rhs->getOpPtrVector().push_back(
new OpSpatialPhysical(stretchTensor, dataAtPts, alphaU));
}
fe_rhs->getOpPtrVector().push_back(
new OpSpatialConsistencyP(piolaStress, dataAtPts));
fe_rhs->getOpPtrVector().push_back(
new OpSpatialConsistencyBubble(bubbleField, dataAtPts));
fe_rhs->getOpPtrVector().push_back(
new OpSpatialConsistencyDivTerm(piolaStress, dataAtPts));
using BodyNaturalBC =
Assembly<PETSC>::LinearForm<
GAUSS>;
BodyNaturalBC::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, 3>;
CHKERR BodyNaturalBC::AddFluxToPipeline<OpBodyForce>::add(
fe_rhs->getOpPtrVector(), mField, spatialL2Disp, {time_scale},
"BODY_FORCE", Sev::inform);
}
CHKERR setBaseVolumeElementOps(tag,
true,
true, fe_lhs);
if (add_elastic) {
const bool symmetric_system =
if(noStretch) {
fe_lhs->getOpPtrVector().push_back(
new OpSpatialConsistency_dP_dP(piolaStress, piolaStress, dataAtPts));
fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dBubble_dP(
bubbleField, piolaStress, dataAtPts));
fe_lhs->getOpPtrVector().push_back(
new OpSpatialConsistency_dBubble_dBubble(bubbleField, bubbleField,
dataAtPts));
} else {
fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_du(
stretchTensor, stretchTensor, dataAtPts, alphaU));
fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dP(
stretchTensor, piolaStress, dataAtPts, true));
fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dBubble(
stretchTensor, bubbleField, dataAtPts, true));
if (!symmetric_system) {
fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_domega(
stretchTensor, rotAxis, dataAtPts, false));
}
}
fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dP(
spatialL2Disp, piolaStress, dataAtPts, true));
fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dw(
spatialL2Disp, spatialL2Disp, dataAtPts, alphaW, alphaRho));
fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dP_domega(
piolaStress, rotAxis, dataAtPts, symmetric_system));
fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dBubble_domega(
bubbleField, rotAxis, dataAtPts, symmetric_system));
if (!symmetric_system) {
fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dP(
rotAxis, piolaStress, dataAtPts, false));
fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dBubble(
rotAxis, bubbleField, dataAtPts, false));
fe_lhs->getOpPtrVector().push_back(
new OpSpatialRotation_domega_domega(rotAxis, rotAxis, dataAtPts));
}
Assembly<BLOCK_PRECONDITIONER_SCHUR>::BiLinearForm<
GAUSS>;
using OpMassStab = B::OpMassCache<1, SPACE_DIM>;
if (precEpsOmega > std::numeric_limits<double>::epsilon()) {
auto cache = boost::make_shared<MatrixDouble>();
fe_lhs->getOpPtrVector().push_back(
new OpMassStab(cache_mats, rotAxis, rotAxis, precEpsOmega));
}
if (std::abs(alphaRho + alphaW) <
std::numeric_limits<double>::epsilon()) {
auto cache = boost::make_shared<MatrixDouble>();
fe_lhs->getOpPtrVector().push_back(
new OpMassStab(cache_mats, spatialL2Disp, spatialL2Disp, precEpsW));
}
}
if (add_material) {
}
}
}
const bool add_elastic, const bool add_material,
boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_rhs,
boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_lhs) {
fe_rhs = boost::make_shared<FaceElementForcesAndSourcesCore>(mField);
fe_lhs = boost::make_shared<FaceElementForcesAndSourcesCore>(mField);
{HDIV});
{HDIV});
if (add_elastic) {
fe_rhs->getOpPtrVector().push_back(
new OpDispBc(piolaStress, dataAtPts, bcSpatialDispVecPtr,
{
boost::make_shared<TimeScale>("disp_history.txt")
}));
fe_rhs->getOpPtrVector().push_back(
new OpRotationBc(piolaStress, dataAtPts, bcSpatialRotationVecPtr));
}
}
boost::shared_ptr<ContactTree> &fe_contact_tree,
boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_face_rhs,
boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_face_lhs
) {
auto get_body_range = [this](auto name, int dim) {
std::map<int, Range> map;
for (auto m_ptr :
(boost::format("%s(.*)") % name).str()
))
) {
CHK_MOAB_THROW(m_ptr->getMeshsetIdEntitiesByDimension(mField.get_moab(),
dim, ents, true),
"by dim");
map[m_ptr->getMeshsetId()] = ents;
}
return map;
};
auto get_map_skin = [this](auto &&map) {
ParallelComm *pcomm =
Skinner skin(&mField.get_moab());
CHKERR skin.find_skin(0,
m.second,
false, skin_faces);
PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, nullptr),
"filter");
m.second.swap(skin_faces);
}
return map;
};
fe_face_rhs = boost::make_shared<BoundaryEle>(mField);
fe_face_lhs = boost::make_shared<BoundaryEle>(mField);
auto rule = [](
int,
int,
int o) {
return -1; };
int levels = 0;
&levels, PETSC_NULL);
auto refine = Tools::refineTriangle(levels);
auto set_rule = [levels, refine](
int order_col, int order_data
) {
auto rule = 2 * order_data;
fe_raw_ptr->gaussPts = Tools::refineTriangleIntegrationPts(rule, refine);
};
fe_face_rhs->getRuleHook = rule;
fe_face_lhs->getRuleHook = rule;
fe_face_rhs->setRuleHook = set_rule;
fe_face_lhs->setRuleHook = set_rule;
fe_face_rhs->getOpPtrVector(), {HDIV});
fe_face_lhs->getOpPtrVector(), {HDIV});
auto add_contact_three = [&]() {
auto tree_moab_ptr = boost::make_shared<moab::Core>();
fe_contact_tree = boost::make_shared<ContactTree>(
mField, tree_moab_ptr, spaceOrder,
get_map_skin(get_body_range("BODY", 3)));
auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
fe_contact_tree->getOpPtrVector(), {HDIV});
fe_contact_tree->getOpPtrVector().push_back(
contactDisp, contact_common_data_ptr->contactDispPtr()));
fe_contact_tree->getOpPtrVector().push_back(
piolaStress, contact_common_data_ptr->contactTractionPtr()));
auto u_h1_ptr = boost::make_shared<MatrixDouble>();
fe_contact_tree->getOpPtrVector().push_back(
fe_contact_tree->getOpPtrVector().push_back(
new OpMoveNode(fe_contact_tree, contact_common_data_ptr, u_h1_ptr));
};
auto add_ops_face_lhs = [&](auto &pip) {
auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
contactDisp, contact_common_data_ptr->contactDispPtr()));
piolaStress, contact_common_data_ptr->contactTractionPtr()));
auto u_h1_ptr = boost::make_shared<MatrixDouble>();
pip.push_back(new OpTreeSearch(
fe_contact_tree, contact_common_data_ptr, u_h1_ptr,
auto contact_sfd_map_range_ptr = boost::make_shared<std::map<int, Range>>(
get_body_range("CONTACT_SDF", 2));
contactDisp, contactDisp, contact_common_data_ptr, fe_contact_tree,
contact_sfd_map_range_ptr));
contactDisp, piolaStress, contact_common_data_ptr, fe_contact_tree,
contact_sfd_map_range_ptr));
piolaStress, contactDisp, contact_common_data_ptr, fe_contact_tree));
using OpMassStab = B::OpMass<1, SPACE_DIM>;
if (precEpsContactDisp > std::numeric_limits<double>::epsilon()) {
pip.push_back(new OpMassStab(
contactDisp, contactDisp,
[this](double, double, double) { return precEpsContactDisp; }));
}
}
};
auto add_ops_face_rhs = [&](auto &pip) {
auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
piolaStress, contact_common_data_ptr->contactTractionPtr()));
contactDisp, contact_common_data_ptr->contactDispPtr()));
auto u_h1_ptr = boost::make_shared<MatrixDouble>();
pip.push_back(new OpTreeSearch(
fe_contact_tree, contact_common_data_ptr, u_h1_ptr,
auto contact_sfd_map_range_ptr = boost::make_shared<std::map<int, Range>>(
get_body_range("CONTACT_SDF", 2));
contactDisp, contact_common_data_ptr, fe_contact_tree,
contact_sfd_map_range_ptr));
piolaStress, contact_common_data_ptr, fe_contact_tree));
};
CHKERR add_ops_face_lhs(fe_face_lhs->getOpPtrVector());
CHKERR add_ops_face_rhs(fe_face_rhs->getOpPtrVector());
}
CHKERR setVolumeElementOps(tag,
true,
false, elasticFeRhs, elasticFeLhs);
CHKERR setFaceElementOps(
true,
false, elasticBcRhs, elasticBcLhs);
auto adj_cache =
boost::make_shared<ForcesAndSourcesCore::UserDataOperator::AdjCache>();
auto get_op_contact_bc = [&]() {
mField, contactElement,
SPACE_DIM - 1, Sev::noisy, adj_cache);
return op_loop_side;
};
auto op_contact_bc = get_op_contact_bc();
elasticFeLhs->getOpPtrVector().push_back(op_contact_bc);
contactTreeRhs, contactRhs, op_contact_bc->getSideFEPtr()
);
}
boost::shared_ptr<FEMethod> null;
boost::shared_ptr<FeTractionBc> spatial_traction_bc(
new FeTractionBc(mField, piolaStress, bcSpatialTraction));
if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
null);
null);
null);
spatial_traction_bc);
null);
null);
} else {
null);
null);
null);
spatial_traction_bc);
null);
null);
}
auto set_set_essential = [&]() {
auto fe_rhs_ptr = boost::make_shared<FEMethod>();
auto fe_lhs_ptr = boost::make_shared<FEMethod>();
auto get_dofs = [&](auto fe_ptr) {
auto vec_glob_dofs =
boost::shared_ptr<std::vector<int>>(fe_ptr, new std::vector<int>());
const auto field_id = mField.get_field_bit_number(piolaStress);
auto numered_dofs_ptr = prb_ptr->getNumeredRowDofsPtr();
vec_glob_dofs->reserve(numered_dofs_ptr->size());
for (
auto pit = (*bcSpatialFreeTraction)[
d].pair_begin();
pit != (*bcSpatialFreeTraction)[
d].pair_end(); ++pit) {
DofEntity::getLoFieldEntityUId(field_id, pit->first));
DofEntity::getHiFieldEntityUId(field_id, pit->second));
for (; lo != hi; ++lo) {
if ((*lo)->getDofCoeffIdx() ==
d) {
vec_glob_dofs->push_back((*lo)->getPetscGlobalDofIdx());
}
}
}
}
return vec_glob_dofs;
};
fe_rhs_ptr->postProcessHook = [this, fe_rhs_ptr]() {
(*fe_rhs_ptr->vecAssembleSwitch) = false;
CHKERR VecGhostUpdateBegin(fe_rhs_ptr->f, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecGhostUpdateEnd(fe_rhs_ptr->f, ADD_VALUES, SCATTER_REVERSE);
CHKERR VecAssemblyBegin(fe_rhs_ptr->f);
CHKERR VecAssemblyEnd(fe_rhs_ptr->f);
CHKERR VecGetArray(fe_rhs_ptr->f, &
a);
const auto field_id = mField.get_field_bit_number(piolaStress);
auto numered_dofs_ptr = fe_rhs_ptr->problemPtr->getNumeredRowDofsPtr();
for (
auto pit = (*bcSpatialFreeTraction)[
d].pair_begin();
pit != (*bcSpatialFreeTraction)[
d].pair_end(); ++pit) {
DofEntity::getLoFieldEntityUId(field_id, pit->first));
DofEntity::getHiFieldEntityUId(field_id, pit->second));
for (; lo != hi; ++lo) {
if ((*lo)->getDofCoeffIdx() ==
d) {
a[(*lo)->getPetscLocalDofIdx()] = 0.0;
}
}
}
}
CHKERR VecRestoreArray(fe_rhs_ptr->f, &
a);
};
auto glob_dofs = get_dofs(fe_lhs_ptr);
fe_lhs_ptr->postProcessHook = [this, fe_lhs_ptr, glob_dofs]() {
(*fe_lhs_ptr->matAssembleSwitch) = false;
CHKERR MatAssemblyBegin(fe_lhs_ptr->B, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(fe_lhs_ptr->B, MAT_FINAL_ASSEMBLY);
CHKERR MatZeroRowsColumns(fe_lhs_ptr->B, glob_dofs->size(),
&*glob_dofs->begin(), 1.0, PETSC_NULL,
PETSC_NULL);
};
ts_ctx_ptr->getPostProcessIFunction().push_back(fe_rhs_ptr);
ts_ctx_ptr->getPostProcessIJacobian().push_back(fe_lhs_ptr);
};
}
#ifdef PYTHON_SDF
boost::shared_ptr<ContactOps::SDFPython> sdf_python_ptr;
auto file_exists = [](std::string myfile) {
std::ifstream file(myfile.c_str());
if (file) {
return true;
}
return false;
};
if (file_exists("sdf.py")) {
MOFEM_LOG(
"EP", Sev::inform) <<
"sdf.py file found";
sdf_python_ptr = boost::make_shared<ContactOps::SDFPython>();
CHKERR sdf_python_ptr->sdfInit(
"sdf.py");
ContactOps::sdfPythonWeakPtr = sdf_python_ptr;
} else {
MOFEM_LOG(
"EP", Sev::warning) <<
"sdf.py file NOT found";
}
#else
#endif
boost::shared_ptr<TsCtx>
ts_ctx;
CHKERR TSAppendOptionsPrefix(ts,
"elastic_");
CHKERR TSSetDM(ts, dmElastic);
SNES snes;
PetscViewerAndFormat *vf;
CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
PETSC_VIEWER_DEFAULT, &vf);
snes,
(
MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
void *))SNESMonitorFields,
PetscSection section;
CHKERR DMGetSection(dmElastic, §ion);
int num_fields;
CHKERR PetscSectionGetNumFields(section, &num_fields);
for (int ff = 0; ff != num_fields; ff++) {
}
CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
boost::shared_ptr<SetUpSchur> schur_ptr;
}
if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
CHKERR TS2SetSolution(ts, x, xx);
} else {
}
TetPolynomialBase::swichCacheHDivBaseDemkowiczOn(
{elasticFeLhs.get(), elasticFeRhs.get()});
CHKERR TSSolve(ts, PETSC_NULL);
TetPolynomialBase::swichCacheHDivBaseDemkowiczOff(
{elasticFeLhs.get(), elasticFeRhs.get()});
int lin_solver_iterations;
CHKERR SNESGetLinearSolveIterations(snes, &lin_solver_iterations);
<< "Number of linear solver iterations " << lin_solver_iterations;
PetscBool test_cook_flg = PETSC_FALSE;
PETSC_NULL);
if (test_cook_flg) {
constexpr int expected_lin_solver_iterations = 11;
}
}
const std::string file) {
if (!dataAtPts) {
dataAtPts =
boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
}
{HDIV});
auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
auto domain_ops = [&](auto &fe) {
fe.getUserPolynomialBase() =
boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
{HDIV, H1, L2});
piolaStress, dataAtPts->getApproxPAtPts()));
bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
stretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
spatialL2Disp, dataAtPts->getSmallWL2AtPts(), MBTET));
spatialH1Disp, dataAtPts->getSmallWH1AtPts()));
spatialH1Disp, dataAtPts->getSmallWGradH1AtPts()));
fe.getOpPtrVector().push_back(
new OpCalculateRotationAndSpatialGradient(dataAtPts));
fe.getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
tag, true, false, dataAtPts, physicalEquations));
spatialL2Disp, contact_common_data_ptr->contactDispPtr()));
};
CHKERR domain_ops(*(op_loop_side->getSideFEPtr()));
post_proc.getOpPtrVector().push_back(op_loop_side);
contactDisp, dataAtPts->getContactL2AtPts()));
post_proc.getOpPtrVector().push_back(new OpPostProcDataStructure(
post_proc.getPostProcMesh(), post_proc.getMapGaussPts(), dataAtPts));
auto u_h1_ptr = boost::make_shared<MatrixDouble>();
auto lambda_h1_ptr = boost::make_shared<MatrixDouble>();
post_proc.getOpPtrVector().push_back(
post_proc.getOpPtrVector().push_back(
piolaStress, contact_common_data_ptr->contactTractionPtr()));
post_proc.getOpPtrVector().push_back(new OpTreeSearch(
contactTreeRhs, contact_common_data_ptr, u_h1_ptr,
&post_proc.getMapGaussPts()));
CHKERR post_proc.writeFile(file.c_str());
}
auto post_proc_norm_fe =
boost::make_shared<VolumeElementForcesAndSourcesCore>(mField);
auto post_proc_norm_rule_hook = [](
int,
int,
int p) ->
int {
return 2 * (p);
};
post_proc_norm_fe->getRuleHook = post_proc_norm_rule_hook;
post_proc_norm_fe->getUserPolynomialBase() =
boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
post_proc_norm_fe->getOpPtrVector(), {L2, H1, HDIV});
enum NORMS { U_NORM_L2 = 0, U_NORM_H1, PIOLA_NORM, U_ERROR_L2, LAST_NORM };
auto norms_vec =
CHKERR VecZeroEntries(norms_vec);
auto u_l2_ptr = boost::make_shared<MatrixDouble>();
auto u_h1_ptr = boost::make_shared<MatrixDouble>();
post_proc_norm_fe->getOpPtrVector().push_back(
post_proc_norm_fe->getOpPtrVector().push_back(
post_proc_norm_fe->getOpPtrVector().push_back(
post_proc_norm_fe->getOpPtrVector().push_back(
post_proc_norm_fe->getOpPtrVector().push_back(
u_h1_ptr));
auto piola_ptr = boost::make_shared<MatrixDouble>();
post_proc_norm_fe->getOpPtrVector().push_back(
post_proc_norm_fe->getOpPtrVector().push_back(
TetPolynomialBase::swichCacheHDivBaseDemkowiczOn({post_proc_norm_fe.get()});
CHKERR mField.loop_finite_elements(
"ELASTIC_PROBLEM", elementVolumeName,
*post_proc_norm_fe);
TetPolynomialBase::swichCacheHDivBaseDemkowiczOff({post_proc_norm_fe.get()});
CHKERR VecAssemblyBegin(norms_vec);
CHKERR VecAssemblyEnd(norms_vec);
const double *norms;
CHKERR VecGetArrayRead(norms_vec, &norms);
MOFEM_LOG(
"EP", Sev::inform) <<
"norm_u: " << std::sqrt(norms[U_NORM_L2]);
MOFEM_LOG(
"EP", Sev::inform) <<
"norm_u_h1: " << std::sqrt(norms[U_NORM_H1]);
<< "norm_error_u_l2: " << std::sqrt(norms[U_ERROR_L2]);
<< "norm_piola: " << std::sqrt(norms[PIOLA_NORM]);
CHKERR VecRestoreArrayRead(norms_vec, &norms);
}
auto bc_mng = mField.getInterface<
BcManager>();
"", piolaStress, false, false);
bcSpatialDispVecPtr = boost::make_shared<BcDispVec>();
for (auto bc : bc_mng->getBcMapByBlockName()) {
if (auto disp_bc = bc.second->dispBcPtr) {
std::vector<double> block_attributes(6, 0.);
if (disp_bc->data.flag1 == 1) {
block_attributes[0] = disp_bc->data.value1;
block_attributes[3] = 1;
}
if (disp_bc->data.flag2 == 1) {
block_attributes[1] = disp_bc->data.value2;
block_attributes[4] = 1;
}
if (disp_bc->data.flag3 == 1) {
block_attributes[2] = disp_bc->data.value3;
block_attributes[5] = 1;
}
auto faces = bc.second->bcEnts.subset_by_dimension(2);
bcSpatialDispVecPtr->emplace_back(bc.first, block_attributes, faces);
}
}
CHKERR getBc(bcSpatialDispVecPtr,
"SPATIAL_DISP_BC", 6);
}
auto bc_mng = mField.getInterface<
BcManager>();
false, false);
bcSpatialTraction = boost::make_shared<TractionBcVec>();
for (auto bc : bc_mng->getBcMapByBlockName()) {
if (auto force_bc = bc.second->forceBcPtr) {
std::vector<double> block_attributes(6, 0.);
block_attributes[0] = -force_bc->data.value3 * force_bc->data.value1;
block_attributes[3] = 1;
block_attributes[1] = -force_bc->data.value4 * force_bc->data.value1;
block_attributes[4] = 1;
block_attributes[2] = -force_bc->data.value5 * force_bc->data.value1;
block_attributes[5] = 1;
auto faces = bc.second->bcEnts.subset_by_dimension(2);
bcSpatialTraction->emplace_back(bc.first, block_attributes, faces);
}
}
CHKERR getBc(bcSpatialTraction,
"SPATIAL_TRACTION_BC", 6);
}
}