#ifndef EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION 3
#endif
OpSetContravariantPiolaTransformOnEdge2D;
};
OpHOSetContravariantPiolaTransformOnFace3D;
};
using OpMixDivULhs = FormsIntegrators<DomainEleOp>::Assembly<
GAUSS>::OpMixDivTimesU<3, SPACE_DIM, SPACE_DIM>;
using OpSpringLhs = FormsIntegrators<BoundaryEleOp>::Assembly<
using OpSpringRhs = FormsIntegrators<BoundaryEleOp>::Assembly<
GAUSS>::OpSource<1, SPACE_DIM>;
GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
GAUSS>::OpGradTensorGrad<1, SPACE_DIM, SPACE_DIM, 1>;
constexpr
double rho = 0;
constexpr
double cn = 0.01;
private:
boost::shared_ptr<ContactOps::CommonData> commonDataPtr;
boost::shared_ptr<PostProcEle> postProcFe;
std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
template <int DIM> Range getEntsOnMeshSkin();
};
}
Simple *
simple = mField.getInterface<Simple>();
enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz"};
PetscInt choice_base_value = AINSWORTH;
LASBASETOPT, &choice_base_value, PETSC_NULL);
switch (choice_base_value) {
case AINSWORTH:
<< "Set AINSWORTH_LEGENDRE_BASE for displacements";
break;
case DEMKOWICZ:
<< "Set DEMKOWICZ_JACOBI_BASE for displacents";
break;
default:
break;
}
auto skin_edges = getEntsOnMeshSkin<SPACE_DIM>();
Range boundary_ents;
ParallelComm *pcomm =
if (pcomm == NULL) {
"Communicator not created");
}
CHKERR pcomm->filter_pstatus(skin_edges, PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &boundary_ents);
}
auto set_matrial_stiffness = [&]() {
: 1;
auto t_D =
getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*commonDataPtr->mDPtr);
};
commonDataPtr = boost::make_shared<ContactOps::CommonData>();
commonDataPtr->mGradPtr = boost::make_shared<MatrixDouble>();
commonDataPtr->mStrainPtr = boost::make_shared<MatrixDouble>();
commonDataPtr->mStressPtr = boost::make_shared<MatrixDouble>();
commonDataPtr->contactStressPtr = boost::make_shared<MatrixDouble>();
commonDataPtr->contactStressDivergencePtr =
boost::make_shared<MatrixDouble>();
commonDataPtr->contactTractionPtr = boost::make_shared<MatrixDouble>();
commonDataPtr->contactDispPtr = boost::make_shared<MatrixDouble>();
commonDataPtr->curlContactStressPtr = boost::make_shared<MatrixDouble>();
commonDataPtr->mDPtr = boost::make_shared<MatrixDouble>();
commonDataPtr->mDPtr->resize(size_symm * size_symm, 1);
CHKERR set_matrial_stiffness();
}
auto bc_mng = mField.getInterface<BcManager>();
auto simple = mField.getInterface<Simple>();
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"NO_CONTACT", "SIGMA", 0, 3);
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_X",
"U", 0, 0);
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_Y",
"U", 1, 1);
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_Z",
"U", 2, 2);
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"REMOVE_ALL", "U", 0, 3);
CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"FIX_X",
"U",
0, 0);
CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"FIX_Y",
"U",
1, 1);
CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"FIX_Z",
"U",
2, 2);
CHKERR bc_mng->pushMarkDOFsOnEntities(
simple->getProblemName(),
"FIX_ALL",
"U", 0, 3);
boundaryMarker = bc_mng->getMergedBlocksMarker(vector<string>{"FIX_"});
}
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
auto bc_mng = mField.getInterface<BcManager>();
auto add_domain_base_ops = [&](auto &pipeline) {
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
pipeline.push_back(new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline.push_back(new OpSetInvJacH1ForFace(inv_jac_ptr));
pipeline.push_back(new OpMakeHdivFromHcurl());
pipeline.push_back(new OpSetHOWeightsOnFace());
} else {
pipeline.push_back(new OpCalculateHOJacVolume(jac_ptr));
pipeline.push_back(new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline.push_back(
new OpSetHOContravariantPiolaTransform(
HDIV, det_ptr, jac_ptr));
pipeline.push_back(
new OpSetHOInvJacVectorBase(
HDIV, inv_jac_ptr));
pipeline.push_back(
new OpSetHOInvJacToScalarBases(
H1, inv_jac_ptr));
pipeline.push_back(new OpSetHOWeights(det_ptr));
}
};
auto henky_common_data_ptr = boost::make_shared<HenckyOps::CommonData>();
henky_common_data_ptr->matGradPtr = commonDataPtr->mGradPtr;
henky_common_data_ptr->matDPtr = commonDataPtr->mDPtr;
auto add_domain_ops_lhs = [&](auto &pipeline) {
pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
pipeline_mng->getOpDomainLhsPipeline().push_back(
"U", commonDataPtr->mGradPtr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpCalculateEigenVals<SPACE_DIM>("U", henky_common_data_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpCalculateLogC<SPACE_DIM>("U", henky_common_data_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpCalculateLogC_dC<SPACE_DIM>("U", henky_common_data_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpCalculateHenckyStress<SPACE_DIM>("U", henky_common_data_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpCalculatePiolaStress<SPACE_DIM>("U", henky_common_data_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpHenckyTangent<SPACE_DIM>("U", henky_common_data_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpKPiola(
"U",
"U", henky_common_data_ptr->getMatTangent()));
} else {
pipeline.push_back(
new OpKCauchy(
"U",
"U", commonDataPtr->mDPtr));
}
auto get_rho = [this](const double, const double, const double) {
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto &fe_domain_lhs = pipeline_mng->getDomainLhsFE();
return rho * fe_domain_lhs->ts_aa;
};
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpMass(
"U",
"U", get_rho));
}
auto unity = []() { return 1; };
pipeline.push_back(
new OpMixDivULhs(
"SIGMA",
"U", unity,
true));
pipeline.push_back(new OpUnSetBc("U"));
};
auto add_domain_ops_rhs = [&](auto &pipeline) {
auto get_body_force = [this](const double, const double, const double) {
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto fe_domain_rhs = pipeline_mng->getDomainRhsFE();
const auto time = fe_domain_rhs->ts_t;
t_source(1) = 1.0 * time;
return t_source;
};
pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
pipeline.push_back(
new OpBodyForce(
"U", get_body_force));
"U", commonDataPtr->mGradPtr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateEigenVals<SPACE_DIM>("U", henky_common_data_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateLogC<SPACE_DIM>("U", henky_common_data_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateLogC_dC<SPACE_DIM>("U", henky_common_data_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateHenckyStress<SPACE_DIM>("U", henky_common_data_ptr));
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculatePiolaStress<SPACE_DIM>("U", henky_common_data_ptr));
"U", henky_common_data_ptr->getMatFirstPiolaStress()));
} else {
pipeline.push_back(new OpSymmetrizeTensor<SPACE_DIM>(
"U", commonDataPtr->mGradPtr, commonDataPtr->mStrainPtr));
pipeline.push_back(new OpTensorTimesSymmetricTensor<SPACE_DIM, SPACE_DIM>(
"U", commonDataPtr->mStrainPtr, commonDataPtr->mStressPtr,
commonDataPtr->mDPtr));
pipeline.push_back(
}
pipeline.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
"U", commonDataPtr->contactDispPtr));
pipeline.push_back(new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(
"SIGMA", commonDataPtr->contactStressPtr));
pipeline.push_back(
new OpCalculateHVecTensorDivergence<SPACE_DIM, SPACE_DIM>(
"SIGMA", commonDataPtr->contactStressDivergencePtr));
pipeline.push_back(
[](double, double, double) { return 1; }));
pipeline.push_back(
"U", commonDataPtr->contactStressDivergencePtr));
pipeline.push_back(
auto mat_acceleration = boost::make_shared<MatrixDouble>();
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpCalculateVectorFieldValuesDotDot<SPACE_DIM>("U",
mat_acceleration));
"U", mat_acceleration, [](
double,
double,
double) {
return rho; }));
}
pipeline.push_back(new OpUnSetBc("U"));
};
auto add_boundary_base_ops = [&](auto &pipeline) {
pipeline.push_back(new OpSetHOWeightsOnFace());
pipeline.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
"U", commonDataPtr->contactDispPtr));
pipeline.push_back(new OpCalculateHVecTensorTrace<SPACE_DIM, BoundaryEleOp>(
"SIGMA", commonDataPtr->contactTractionPtr));
};
auto add_boundary_ops_lhs = [&](auto &pipeline) {
auto &bc_map = bc_mng->getBcMapByBlockName();
for (auto bc : bc_map) {
if (bc_mng->checkBlock(bc, "FIX_")) {
<< "Set boundary matrix for " << bc.first;
pipeline.push_back(
new OpSetBc("U", false, bc.second->getBcMarkersPtr()));
"U", "U", [](double, double, double) { return 1.; },
bc.second->getBcEntsPtr()));
pipeline.push_back(new OpUnSetBc("U"));
}
}
pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
pipeline.push_back(
new OpConstrainBoundaryLhs_dU("SIGMA", "U", commonDataPtr));
pipeline.push_back(
new OpConstrainBoundaryLhs_dTraction("SIGMA", "SIGMA", commonDataPtr));
"U", "U",
));
if (boundaryMarker)
pipeline.push_back(new OpUnSetBc("U"));
};
auto time_scaled = [&](double, double, double) {
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto &fe_domain_rhs = pipeline_mng->getDomainRhsFE();
return -fe_domain_rhs->ts_t;
};
auto add_boundary_ops_rhs = [&](auto &pipeline) {
for (auto &bc : mField.getInterface<BcManager>()->getBcMapByBlockName()) {
if (bc_mng->checkBlock(bc, "FIX_")) {
<< "Set boundary residual for " << bc.first;
pipeline.push_back(
new OpSetBc("U", false, bc.second->getBcMarkersPtr()));
auto attr_vec = boost::make_shared<MatrixDouble>(
SPACE_DIM, 1);
attr_vec->clear();
if (bc.second->bcAttributes.size() !=
SPACE_DIM)
"Wrong size of boundary attributes vector. Set right block "
"size attributes. Size of attributes %d",
bc.second->bcAttributes.size());
std::copy(&bc.second->bcAttributes[0],
attr_vec->data().begin());
bc.second->getBcEntsPtr()));
"U", commonDataPtr->contactDispPtr,
[](double, double, double) { return 1.; },
bc.second->getBcEntsPtr()));
pipeline.push_back(new OpUnSetBc("U"));
}
}
pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
pipeline.push_back(new OpConstrainBoundaryRhs("SIGMA", commonDataPtr));
"U", commonDataPtr->contactDispPtr,
[this](double, double, double) { return spring_stiffness; }));
pipeline.push_back(new OpUnSetBc("U"));
};
add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
add_domain_ops_lhs(pipeline_mng->getOpDomainLhsPipeline());
add_domain_ops_rhs(pipeline_mng->getOpDomainRhsPipeline());
add_boundary_base_ops(pipeline_mng->getOpBoundaryLhsPipeline());
add_boundary_base_ops(pipeline_mng->getOpBoundaryRhsPipeline());
CHKERR add_boundary_ops_lhs(pipeline_mng->getOpBoundaryLhsPipeline());
CHKERR add_boundary_ops_rhs(pipeline_mng->getOpBoundaryRhsPipeline());
};
CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule_vol);
CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule_vol);
};
CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule_boundary);
CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule_boundary);
}
Simple *
simple = mField.getInterface<Simple>();
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
ISManager *is_manager = mField.getInterface<ISManager>();
auto set_section_monitor = [&](auto solver) {
SNES snes;
CHKERR TSGetSNES(solver, &snes);
PetscViewerAndFormat *vf;
CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
PETSC_VIEWER_DEFAULT, &vf);
snes,
(
MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
void *))SNESMonitorFields,
};
auto scatter_create = [&](
auto D,
auto coeff) {
SmartPetscObj<IS> is;
CHKERR is_manager->isCreateProblemFieldAndRank(
simple->getProblemName(),
ROW,
"U", coeff, coeff, is);
int loc_size;
CHKERR ISGetLocalSize(is, &loc_size);
CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &
v);
VecScatter scatter;
CHKERR VecScatterCreate(
D, is,
v, PETSC_NULL, &scatter);
return std::make_tuple(SmartPetscObj<Vec>(
v),
SmartPetscObj<VecScatter>(scatter));
};
auto set_time_monitor = [&](auto dm, auto solver) {
boost::shared_ptr<Monitor> monitor_ptr(
new Monitor(dm, commonDataPtr, uXScatter, uYScatter, uZScatter));
boost::shared_ptr<ForcesAndSourcesCore> null;
monitor_ptr, null, null);
};
uXScatter = scatter_create(
D, 0);
uYScatter = scatter_create(
D, 1);
uZScatter = scatter_create(
D, 2);
auto solver = pipeline_mng->createTSIM();
CHKERR set_section_monitor(solver);
CHKERR set_time_monitor(dm, solver);
CHKERR TSSetFromOptions(solver);
} else {
auto solver = pipeline_mng->createTSIM2();
CHKERR set_section_monitor(solver);
CHKERR set_time_monitor(dm, solver);
CHKERR TS2SetSolution(solver,
D, DD);
CHKERR TSSetFromOptions(solver);
}
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
}
Range body_ents;
CHKERR mField.get_moab().get_entities_by_dimension(0, DIM, body_ents);
Skinner skin(&mField.get_moab());
Range skin_ents;
CHKERR skin.find_skin(0, body_ents,
false, skin_ents);
return skin_ents;
};
static char help[] =
"...\n\n";
int main(
int argc,
char *argv[]) {
auto core_log = logging::core::get();
core_log->add_sink(
try {
DMType dm_name = "DMMOFEM";
}
}
ForcesAndSourcesCore::UserDataOperator UserDataOperator
void simple(double P1[], double P2[], double P3[], double c[], const int N)
int main(int argc, char *argv[])
EntitiesFieldData::EntData EntData
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Kronecker Delta class symmetric.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
constexpr double poisson_ratio
constexpr double shear_modulus_G
constexpr double bulk_modulus_K
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, SPACE_DIM > OpBodyForce
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpInertiaForce
constexpr bool is_quasi_static
constexpr double young_modulus
#define CATCH_ERRORS
Catch errors.
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
FieldSpace
approximation spaces
@ HCURL
field with continuous tangents
@ HDIV
field with continuous normal traction
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
auto smartCreateDMVector
Get smart vector from DM.
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
FTensor::Index< 'i', SPACE_DIM > i
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
FormsIntegrators< BoundaryEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMass< 1, 3 > OpSpringLhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 3, 3 > OpMixUTimesLambdaRhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMixDivTimesVec< 3 > OpMixDivULhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpMixVecTimesDivLambda< 3 > OpMixUTimesDivLambdaRhs
FormsIntegrators< BoundaryEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpSpringRhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpMixDivTimesU< 3, 3, 3 > OpMixDivURhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMixTensorTimesGrad< 3 > OpLambdaGraULhs
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
OpSetInvJacHcurlFaceImpl< 2 > OpSetInvJacHcurlFace
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
OpSetContravariantPiolaTransformOnFace2DImpl< 2 > OpSetContravariantPiolaTransformOnFace2D
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
DeprecatedCoreInterface Interface
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
const double D
diffusivity
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
#define EXECUTABLE_DIMENSION
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 0 > OpBoundaryVec
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Body force]
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpBoundaryInternal
PetscBool is_large_strains
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
constexpr EntityType boundary_ent
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Range getEntsOnMeshSkin()
[Check]
MoFEMErrorCode tsSolve()
[Push operators to pipeline]
MoFEMErrorCode checkResults()
[Postprocess results]
MoFEMErrorCode createCommonData()
[Create common data]
MoFEMErrorCode OPs()
[Boundary condition]
MoFEMErrorCode runProblem()
[Create common data]
MoFEMErrorCode postProcess()
[Solve]
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode bC()
[Create common data]
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
default operator for TRI element
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
MoFEM::FaceElementForcesAndSourcesCore FaceEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
[Push operators to pipeline]
boost::shared_ptr<MatrixDouble>
mDPtr;
boost::shared_ptr<MatrixDouble>
mGradPtr;
};
boost::shared_ptr<CommonData> common_data_ptr);
private:
};
const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr);
private:
};
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr);
private:
};
template <typename T1, typename T2>
t_normal(1) = 1.;
return t_normal;
}
template <typename T>
return (-0.5 - t_coords(1)) * t_normal(1);
}
template <typename T>
return t_disp(
i) * t_normal(
i);
}
template <typename T>
return -t_traction(
i) * t_normal(
i);
}
inline double sign(
double x) {
if (x == 0)
return 0;
else if (x > 0)
return 1;
else
return -1;
};
inline double w(
const double g,
const double t) {
return g -
cn *
t; }
inline double constrian(
double &&g0,
double &&
g,
double &&
t) {
return (
w(
g - g0,
t) + std::abs(
w(
g - g0,
t))) / 2 + g0;
};
}
return (1 +
sign(
w(
g - g0,
t))) / 2;
}
const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
commonDataPtr(common_data_ptr) {}
const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
auto &nf = AssemblyBoundaryEleOp::locF;
auto t_normal = getFTensor1Normal();
t_normal(
i) /= sqrt(t_normal(
j) * t_normal(
j));
auto t_w = getFTensor0IntegrationWeight();
auto t_disp = getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->contactDispPtr));
auto t_traction =
getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->contactTractionPtr));
auto t_coords = getFTensor1CoordsAtGaussPts();
size_t nb_base_functions = data.getN().size2() / 3;
auto t_base = data.getFTensor1N<3>();
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
auto t_nf = getFTensor1FromPtr<SPACE_DIM>(&nf[0]);
const double alpha = t_w * getMeasure();
auto t_contact_normal =
normal(t_coords, t_disp);
t_P(
i,
j) = t_contact_normal(
i) * t_contact_normal(
j);
gap(t_disp, t_contact_normal),
t_rhs_tangent_traction;
t_rhs_tangent_disp(
i) = t_Q(
i,
j) * t_disp(
j);
t_rhs_tangent_traction(
i) =
cn * t_Q(
i,
j) * t_traction(
j);
size_t bb = 0;
for (; bb != AssemblyBoundaryEleOp::nbRows /
SPACE_DIM; ++bb) {
const double beta =
alpha * (t_base(
i) * t_normal(
i));
t_nf(
i) -= beta * t_rhs_constrains(
i);
t_nf(
i) -= beta * t_rhs_tangent_disp(
i);
t_nf(
i) += beta * t_rhs_tangent_traction(
i);
++t_nf;
++t_base;
}
for (; bb < nb_base_functions; ++bb)
++t_base;
++t_disp;
++t_traction;
++t_coords;
++t_w;
}
}
OpConstrainBoundaryLhs_dU::OpConstrainBoundaryLhs_dU(
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr)
commonDataPtr(common_data_ptr) {
sYmm = false;
}
const size_t nb_gauss_pts = getGaussPts().size2();
auto &locMat = AssemblyBoundaryEleOp::locMat;
auto t_normal = getFTensor1Normal();
t_normal(
i) /= sqrt(t_normal(
j) * t_normal(
j));
auto t_disp = getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->contactDispPtr));
auto t_traction =
getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->contactTractionPtr));
auto t_coords = getFTensor1CoordsAtGaussPts();
auto t_w = getFTensor0IntegrationWeight();
auto t_row_base = row_data.getFTensor1N<3>();
size_t nb_face_functions = row_data.getN().size2() / 3;
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
const double alpha = t_w * getMeasure();
auto t_contact_normal =
normal(t_coords, t_disp);
t_P(
i,
j) = t_contact_normal(
i) * t_contact_normal(
j);
gap0(t_coords, t_contact_normal),
gap(t_disp, t_contact_normal),
size_t rr = 0;
for (; rr != AssemblyBoundaryEleOp::nbRows /
SPACE_DIM; ++rr) {
auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
const double row_base = t_row_base(
i) * t_normal(
i);
auto t_col_base = col_data.getFTensor0N(gg, 0);
for (
size_t cc = 0; cc != AssemblyBoundaryEleOp::nbCols /
SPACE_DIM;
++cc) {
const double beta =
alpha * row_base * t_col_base;
t_mat(
i,
j) -= (beta * diff_constrain) * t_P(
i,
j);
t_mat(
i,
j) -= beta * t_Q(
i,
j);
++t_col_base;
++t_mat;
}
++t_row_base;
}
for (; rr < nb_face_functions; ++rr)
++t_row_base;
++t_disp;
++t_traction;
++t_coords;
++t_w;
}
}
OpConstrainBoundaryLhs_dTraction::OpConstrainBoundaryLhs_dTraction(
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr)
commonDataPtr(common_data_ptr) {
sYmm = false;
}
const size_t nb_gauss_pts = getGaussPts().size2();
auto &locMat = AssemblyBoundaryEleOp::locMat;
auto t_normal = getFTensor1Normal();
t_normal(
i) /= sqrt(t_normal(
j) * t_normal(
j));
auto t_disp = getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->contactDispPtr));
auto t_traction =
getFTensor1FromMat<SPACE_DIM>(*(commonDataPtr->contactTractionPtr));
auto t_coords = getFTensor1CoordsAtGaussPts();
auto t_w = getFTensor0IntegrationWeight();
auto t_row_base = row_data.getFTensor1N<3>();
size_t nb_face_functions = row_data.getN().size2() / 3;
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
const double alpha = t_w * getMeasure();
auto t_contact_normal =
normal(t_coords, t_disp);
t_P(
i,
j) = t_contact_normal(
i) * t_contact_normal(
j);
gap0(t_coords, t_contact_normal),
gap(t_disp, t_contact_normal),
size_t rr = 0;
for (; rr != AssemblyBoundaryEleOp::nbRows /
SPACE_DIM; ++rr) {
auto t_mat = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
const double row_base = t_row_base(
i) * t_normal(
i);
auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
for (
size_t cc = 0; cc != AssemblyBoundaryEleOp::nbCols /
SPACE_DIM;
++cc) {
const double col_base = t_col_base(
i) * t_normal(
i);
const double beta =
alpha * row_base * col_base;
t_mat(
i,
j) += (beta * diff_traction) * t_P(
i,
j);
t_mat(
i,
j) += beta *
cn * t_Q(
i,
j);
++t_col_base;
++t_mat;
}
++t_row_base;
}
for (; rr < nb_face_functions; ++rr)
++t_row_base;
++t_disp;
++t_traction;
++t_coords;
++t_w;
}
}
};
Tensor2_Expr< Kronecker_Delta< T >, T, Dim0, Dim1, i, j > kronecker_delta(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
Rank 2.
constexpr double t
plate stiffness