#ifndef EXECUTABLE_DIMENSION
#define EXECUTABLE_DIMENSION 3
#endif
#ifndef SCHUR_ASSEMBLE
#define SCHUR_ASSEMBLE 0
#endif
#include <GenericElementInterface.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;
#endif
IntegrationType::GAUSS;
};
};
IT>::OpBaseTimesVector<1, SPACE_DIM, 1>;
};
#ifdef WITH_MODULE_MFRONT_INTERFACE
#include <MFrontMoFEMInterface.hpp>
#endif
private:
#ifdef PYTHON_SDF
boost::shared_ptr<SDFPython> sdfPythonPtr;
#endif
};
};
};
}
PETSC_NULL);
PETSC_NULL);
PETSC_NULL);
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 displacements";
break;
default:
break;
}
CHKERR skin.find_skin(0, body_ents,
false, skin_ents);
return skin_ents;
};
auto filter_blocks = [&](auto skin) {
bool is_contact_block = false;
(boost::format("%s(.*)") % "CONTACT").str()
))
) {
is_contact_block =
true;
<<
"Find contact block set: " <<
m->getName();
auto meshset =
m->getMeshset();
Range contact_meshset_range;
meshset,
SPACE_DIM - 1, contact_meshset_range,
true);
contact_meshset_range);
contact_range.merge(contact_meshset_range);
}
if (is_contact_block) {
<< "Nb entities in contact surface: " << contact_range.size();
skin = intersect(skin, contact_range);
}
return skin;
};
ParallelComm *pcomm =
CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &boundary_ents);
return boundary_ents;
};
moab::Interface::UNION);
}
auto project_ho_geometry = [&]() {
};
PetscBool project_geometry = PETSC_TRUE;
&project_geometry, PETSC_NULL);
if (project_geometry) {
}
}
PetscBool use_mfront = PETSC_FALSE;
PETSC_NULL);
PETSC_NULL);
PETSC_NULL);
PETSC_NULL);
PETSC_NULL);
} else {
MOFEM_LOG(
"CONTACT", Sev::inform) <<
"Using MFront for material model";
}
PetscBool use_scale = PETSC_FALSE;
PETSC_NULL);
if (use_scale) {
}
#ifdef PYTHON_SDF
char sdf_file_name[255];
sdf_file_name, 255, PETSC_NULL);
sdfPythonPtr = boost::make_shared<SDFPython>();
CHKERR sdfPythonPtr->sdfInit(sdf_file_name);
sdfPythonWeakPtr = sdfPythonPtr;
#endif
"Use executable contact_2d with axisymmetric model");
} else {
if (!use_mfront) {
"Axisymmetric model is only available with MFront (set "
"use_mfront to 1)");
} else {
MOFEM_LOG(
"CONTACT", Sev::inform) <<
"Using axisymmetric model";
}
}
} else {
MOFEM_LOG(
"CONTACT", Sev::inform) <<
"Using plane strain model";
}
}
if (use_mfront) {
#ifndef WITH_MODULE_MFRONT_INTERFACE
SETERRQ(
"MFrontInterface module was not found while use_mfront was set to 1");
#else
boost::make_shared<MFrontMoFEMInterface<TRIDIMENSIONAL>>(
boost::make_shared<MFrontMoFEMInterface<AXISYMMETRICAL>>(
} else {
}
}
#endif
}
if (use_mfront) {
}
}
for (
auto f : {
"U",
"SIGMA"}) {
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
}
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"FIX_X",
"SIGMA", 0, 0, false, true);
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"FIX_Y",
"SIGMA", 1, 1, false, true);
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"FIX_Z",
"SIGMA", 2, 2, false, true);
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"FIX_ALL",
"SIGMA", 0, 3, false, true);
CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"NO_CONTACT",
"SIGMA", 0, 3,
false,
true);
simple->getProblemName(),
"U");
}
auto time_scale = boost::make_shared<ScaledTimeScale>();
auto body_force_time_scale =
boost::make_shared<ScaledTimeScale>("body_force_hist.txt");
auto integration_rule_vol = [](int, int,
int approx_order) {
};
auto integration_rule_boundary = [](int, int,
int approx_order) {
};
auto add_domain_base_ops = [&](auto &pip) {
"GEOMETRY");
};
auto hencky_common_data_ptr = boost::make_shared<HenckyOps::CommonData>();
hencky_common_data_ptr->matDPtr = boost::make_shared<MatrixDouble>();
hencky_common_data_ptr->matGradPtr = boost::make_shared<MatrixDouble>();
auto add_domain_ops_lhs = [&](auto &pip) {
auto fe_domain_lhs = pip_mng->getDomainLhsFE();
auto get_inertia_and_mass_damping =
return (
rho *
scale) * fe_domain_lhs->ts_aa +
};
pip.push_back(
new OpMass(
"U",
"U", get_inertia_and_mass_damping));
} else {
auto fe_domain_lhs = pip_mng->getDomainLhsFE();
auto get_inertia_and_mass_damping =
};
pip.push_back(
new OpMass(
"U",
"U", get_inertia_and_mass_damping));
}
} else {
}
};
auto add_domain_ops_rhs = [&](auto &pip) {
pip,
mField,
"U", {body_force_time_scale}, Sev::inform);
auto mat_acceleration = boost::make_shared<MatrixDouble>();
"U", mat_acceleration));
pip.push_back(
new OpInertiaForce(
"U", mat_acceleration, [](
double,
double,
double) {
}));
}
auto mat_velocity = boost::make_shared<MatrixDouble>();
pip.push_back(
pip.push_back(
}));
}
} else {
}
};
auto add_boundary_base_ops = [&](auto &pip) {
"GEOMETRY");
};
auto add_boundary_ops_lhs = [&](auto &pip) {
pip,
mField,
"U", Sev::inform);
auto fe_boundary_lhs = pip_mng->getBoundaryLhsFE();
"U", "U",
[this, fe_boundary_lhs](double, double, double) {
}
));
}
mField, pip,
simple->getDomainFEName(),
"SIGMA",
"U",
"GEOMETRY",
};
auto add_boundary_ops_rhs = [&](auto &pip) {
pip,
mField,
"U", {time_scale}, Sev::inform);
auto u_disp = boost::make_shared<MatrixDouble>();
auto dot_u_disp = boost::make_shared<MatrixDouble>();
pip.push_back(
pip.push_back(
new OpSpringRhs(
"U", u_disp, [
this](
double,
double,
double) {
}));
pip.push_back(
new OpSpringRhs(
"U", dot_u_disp, [
this](
double,
double,
double) {
}));
}
};
CHKERR add_domain_base_ops(pip_mng->getOpDomainLhsPipeline());
CHKERR add_domain_base_ops(pip_mng->getOpDomainRhsPipeline());
CHKERR add_domain_ops_lhs(pip_mng->getOpDomainLhsPipeline());
CHKERR add_domain_ops_rhs(pip_mng->getOpDomainRhsPipeline());
CHKERR add_boundary_base_ops(pip_mng->getOpBoundaryLhsPipeline());
CHKERR add_boundary_base_ops(pip_mng->getOpBoundaryRhsPipeline());
CHKERR add_boundary_ops_lhs(pip_mng->getOpBoundaryLhsPipeline());
CHKERR add_boundary_ops_rhs(pip_mng->getOpBoundaryRhsPipeline());
}
CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule_vol);
CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule_vol);
CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule_boundary);
CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule_boundary);
}
static boost::shared_ptr<SetUpSchur>
protected:
};
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) {
CHKERR is_manager->isCreateProblemFieldAndRank(
simple->getProblemName(),
ROW,
"U", coeff, coeff, is);
int loc_size;
CHKERR ISGetLocalSize(is, &loc_size);
VecScatter scatter;
CHKERR VecScatterCreate(
D, is,
v, PETSC_NULL, &scatter);
};
auto set_time_monitor = [&](auto dm, auto solver) {
boost::shared_ptr<ForcesAndSourcesCore> null;
};
auto set_essential_bc = [&]() {
auto pre_proc_ptr = boost::make_shared<FEMethod>();
auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
auto time_scale = boost::make_shared<TimeScale>();
auto get_bc_hook_rhs = [&]() {
{time_scale}, false);
return hook;
};
pre_proc_ptr->preProcessHook = get_bc_hook_rhs();
auto get_post_proc_hook_rhs = [&]() {
mField, post_proc_rhs_ptr, 1.);
};
auto get_post_proc_hook_lhs = [&]() {
mField, post_proc_lhs_ptr, 1.);
};
post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs();
ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();
ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
};
auto set_schur_pc = [&](auto solver) {
boost::shared_ptr<SetUpSchur> schur_ptr;
if (
AT == AssemblyType::BLOCK_SCHUR) {
}
return schur_ptr;
};
auto solver = pip_mng->createTSIM();
CHKERR TSSetFromOptions(solver);
auto schur_pc_ptr = set_schur_pc(solver);
CHKERR set_section_monitor(solver);
CHKERR set_time_monitor(dm, solver);
} else {
auto solver = pip_mng->createTSIM2();
CHKERR TSSetFromOptions(solver);
auto schur_pc_ptr = set_schur_pc(solver);
CHKERR set_section_monitor(solver);
CHKERR set_time_monitor(dm, solver);
CHKERR TS2SetSolution(solver,
D, DD);
}
}
const double *t_ptr;
double hertz_force;
double fem_force;
double analytical_active_area = 1.0;
double norm = 1e-5;
double tol_force = 1e-3;
double tol_norm = 7.5;
double tol_area = 3e-2;
double fem_active_area = t_ptr[3];
case 1:
hertz_force = 3.927;
fem_force = t_ptr[1];
break;
case 2:
hertz_force = 4.675;
fem_force = t_ptr[1];
break;
case 3:
hertz_force = 3.968;
tol_force = 2e-3;
fem_force = t_ptr[2];
analytical_active_area = M_PI / 4;
tol_area = 0.2;
break;
case 4:
tol_force = 5e-3;
tol_area = 0.2;
case 5:
hertz_force = 15.873;
tol_force = 5e-3;
fem_force = t_ptr[1];
analytical_active_area = M_PI;
break;
case 6:
hertz_force = 0.374;
fem_force = t_ptr[1];
break;
case 7:
hertz_force = 0.5289;
fem_force = t_ptr[2];
break;
default:
}
if (fabs(fem_force - hertz_force) / hertz_force > tol_force) {
"atom test %d failed: Wrong FORCE output: %3.4e != %3.4e",
}
if (norm > tol_norm) {
"atom test %d failed: Wrong NORM output: %3.4e > %3.4e",
}
if (fabs(fem_active_area - analytical_active_area) > tol_area) {
"atom test %d failed: AREA computed %3.4e but should be %3.4e",
atom_test, fem_active_area, analytical_active_area);
}
}
}
static char help[] =
"...\n\n";
int main(
int argc,
char *argv[]) {
#ifdef PYTHON_SDF
Py_Initialize();
np::initialize();
#endif
const char param_file[] = "param_file.petsc";
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "CONTACT"));
LogManager::setLog("CONTACT");
try {
DMType dm_name = "DMMOFEM";
DMType dm_name_mg = "DMMOFEM_MG";
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
}
#ifdef PYTHON_SDF
if (Py_FinalizeEx() < 0) {
exit(120);
}
#endif
return 0;
}
private:
};
SNES snes;
CHKERR TSGetSNES(solver, &snes);
KSP ksp;
CHKERR SNESGetKSP(snes, &ksp);
CHKERR KSPSetFromOptions(ksp);
PC pc;
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
if (is_pcfs) {
MOFEM_LOG(
"CONTACT", Sev::inform) <<
"Setup Schur pc";
"It is expected that Schur matrix is not allocated. This is "
"possible only if PC is set up twice");
}
DM solver_dm;
CHKERR TSGetDM(solver, &solver_dm);
CHKERR DMSetMatType(solver_dm, MATSHELL);
auto swap_assemble = [](TS ts, PetscReal
t, Vec u, Vec u_t, PetscReal
a,
Mat
A, Mat
B,
void *ctx) {
};
CHKERR TSSetIJacobian(solver, A, P, swap_assemble, ts_ctx_ptr.get());
} else {
auto swap_assemble = [](TS ts, PetscReal
t, Vec u, Vec u_t, Vec utt,
PetscReal
a, PetscReal aa, Mat
A, Mat
B,
void *ctx) {
};
CHKERR TSSetI2Jacobian(solver, A, P, swap_assemble, ts_ctx_ptr.get());
}
CHKERR KSPSetOperators(ksp, A, P);
} else {
MOFEM_LOG(
"CONTACT", Sev::inform) <<
"No Schur PC";
}
}
auto create_dm = [&](
const char *name,
const char *
field_name,
auto dm_type) {
auto create_dm_imp = [&]() {
};
create_dm_imp(),
"Error in creating schurDM. It is possible that schurDM is "
"already created");
return dm;
};
schurDM = create_dm(
"SCHUR",
"U",
"DMMOFEM_MG");
blockDM = create_dm(
"BLOCK",
"SIGMA",
"DMMOFEM");
if constexpr (
AT == AssemblyType::BLOCK_SCHUR) {
auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
{{
simple->getDomainFEName(),
{
{"U", "U"}, {"SIGMA", "U"}, {"U", "SIGMA"}, {"SIGMA", "SIGMA"}
}}}
);
{schur_dm, block_dm}, block_mat_data,
{"SIGMA"}, {nullptr}, true
);
};
} else {
"Only BLOCK_SCHUR is implemented");
}
}
double eps_stab = 1e-4;
PETSC_NULL);
using OpMassStab = B::OpMass<3, SPACE_DIM * SPACE_DIM>;
pip->getOpBoundaryLhsPipeline().push_back(
new OpMassStab("SIGMA", "SIGMA",
[eps_stab](double, double, double) { return eps_stab; }));
pip->getOpBoundaryLhsPipeline().push_back(
);
pip->getOpDomainLhsPipeline().push_back(
false)
);
auto pre_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
auto post_proc_schur_lhs_ptr = boost::make_shared<FEMethod>();
pre_proc_schur_lhs_ptr->preProcessHook = [this]() {
MOFEM_LOG(
"CONTACT", Sev::verbose) <<
"Lhs Assemble Begin";
};
post_proc_schur_lhs_ptr->postProcessHook = [this, ao_up,
post_proc_schur_lhs_ptr]() {
MOFEM_LOG(
"CONTACT", Sev::verbose) <<
"Lhs Assemble End";
auto print_mat_norm = [
this](
auto a, std::string prefix) {
double nrm;
CHKERR MatNorm(
a, NORM_FROBENIUS, &nrm);
MOFEM_LOG(
"CONTACT", Sev::noisy) << prefix <<
" norm = " << nrm;
};
CHKERR MatAssemblyBegin(
S, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(
S, MAT_FINAL_ASSEMBLY);
mField, post_proc_schur_lhs_ptr, 1,
S, ao_up)();
#ifndef NDEBUG
#endif
MOFEM_LOG(
"CONTACT", Sev::verbose) <<
"Lhs Assemble Finish";
};
ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_schur_lhs_ptr);
ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_schur_lhs_ptr);
}
CHKERR PCFieldSplitSetIS(pc, NULL, block_is);
CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER,
S);
}
KSP *subksp;
CHKERR PCFieldSplitSchurGetSubKSP(pc, PETSC_NULL, &subksp);
auto get_pc = [](auto ksp) {
PC pc_raw;
CHKERR KSPGetPC(ksp, &pc_raw);
};
auto set_pc_p_mg = [](
auto dm,
auto pc,
auto S) {
PetscBool same = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCMG, &same);
if (same) {
}
};
}
boost::shared_ptr<SetUpSchur>
}
static auto filter_true_skin(MoFEM::Interface &m_field, Range &&skin)
static auto get_skin(MoFEM::Interface &m_field, Range body_ents)
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
void simple(double P1[], double P2[], double P3[], double c[], const int N)
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
[Define dimension]
#define CATCH_ERRORS
Catch errors.
FieldApproximationBase
approximation base
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
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 ...
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
@ MOFEM_ATOM_TEST_INVALID
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpInertiaForce
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
MoFEMErrorCode DMRegister_MGViaApproxOrders(const char sname[])
Register DM for Multi-Grid via approximation orders.
auto createDMVector(DM dm)
Get smart vector from DM.
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
auto createDMMatrix(DM dm)
Get smart matrix from DM.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
const double v
phase velocity of light in medium (cm/ns)
MoFEMErrorCode opFactoryDomainLhs(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, boost::shared_ptr< HenckyOps::CommonData > common_ptr, Sev sev)
MoFEMErrorCode opFactoryDomainRhs(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, boost::shared_ptr< HenckyOps::CommonData > common_ptr, Sev sev)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
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.
PetscErrorCode TsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
Set function evaluating jacobian in TS solver.
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
OpSchurAssembleBase * createOpSchurAssembleEnd(std::vector< std::string > fields_name, std::vector< boost::shared_ptr< Range > > field_ents, SmartPetscObj< AO > ao, SmartPetscObj< Mat > schur, bool sym_schur, bool symm_op)
Construct a new Op Schur Assemble End object.
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
MoFEMErrorCode setSchurA00MatSolvePC(SmartPetscObj< PC > pc)
Set PC for A00 block.
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
boost::shared_ptr< PCMGSetUpViaApproxOrdersCtx > createPCMGSetUpViaApproxOrdersCtx(DM dm, Mat A, bool use_shell_mat)
createPCMGSetUpViaApproxOrdersCtx
auto getDMSubData(DM dm)
Get sub problem data structure.
PetscErrorCode TsSetI2Jacobian(TS ts, PetscReal t, Vec u, Vec u_t, Vec u_tt, PetscReal a, PetscReal aa, Mat A, Mat B, void *ctx)
Calculation Jacobian for second order PDE in time.
MoFEMErrorCode PCMGSetUpViaApproxOrders(PC pc, boost::shared_ptr< PCMGSetUpViaApproxOrdersCtx > ctx, int verb)
Function build MG structure.
boost::shared_ptr< BlockStructure > createBlockMatStructure(DM dm, SchurFEOpsFEandFields schur_fe_op_vec)
Create a Mat Diag Blocks object.
boost::shared_ptr< NestSchurData > createSchurNestedMatrixStruture(std::pair< SmartPetscObj< DM >, SmartPetscObj< DM > > dms, boost::shared_ptr< BlockStructure > block_mat_data_ptr, std::vector< std::string > fields_names, std::vector< boost::shared_ptr< Range > > field_ents, bool add_preconditioner_block)
Get the Schur Nest Mat Array object.
auto createAOMappingIS(IS isapp, IS ispetsc)
Creates an application mapping using two index sets.
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
MoFEMErrorCode DMMoFEMSetNestSchurData(DM dm, boost::shared_ptr< NestSchurData >)
auto createDMNestSchurMat(DM dm)
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
OpSchurAssembleBase * createOpSchurAssembleBegin()
auto createDMBlockMat(DM dm)
double young_modulus
Young modulus.
#define EXECUTABLE_DIMENSION
PetscBool is_quasi_static
double poisson_ratio
Poisson ratio.
int geom_order
Order if fixed.
ElementsAndOps< SPACE_DIM >::SideEle SideEle
constexpr FieldSpace CONTACT_SPACE
constexpr double t
plate stiffness
constexpr auto field_name
static constexpr int approx_order
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
FTensor::Index< 'm', 3 > m
Simple interface for fast problem set-up.
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
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.
Definition of the displacement bc data structure.
Data on single entity (This is passed as argument to DataOperator::doWork)
Class (Function) to enforce essential constrains on the left hand side diagonal.
Class (Function) to enforce essential constrains on the right hand side diagonal.
Class (Function) to enforce essential constrains.
Section manager is used to create indexes and sections.
Interface for managing meshsets containing materials and boundary conditions.
Approximate field values for given petsc vector.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
PipelineManager interface.
Projection of edge entities with one mid-node on hierarchical basis.
Simple interface for fast problem set-up.
intrusive_ptr for managing petsc objects
Force scale operator for reading two columns.
double getScale(const double time)
Get scaling at a given time.
TimeScale(std::string file_name="", bool error_if_file_not_given=false, ScalingFun def_scaling_fun=[](double time) { return time;})
TimeScale constructor.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
[Push operators to pipeline]
virtual MoFEMErrorCode setUp(SmartPetscObj< KSP >)=0
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field)
virtual ~SetUpSchurImpl()=default
MoFEMErrorCode createSubDM()
MoFEMErrorCode setDiagonalPC(PC pc)
SetUpSchurImpl(MoFEM::Interface &m_field)
SmartPetscObj< DM > schurDM
MoFEMErrorCode setUp(SmartPetscObj< KSP >)
SmartPetscObj< DM > blockDM
MoFEMErrorCode setPC(PC pc)
MoFEMErrorCode setOperator()
MoFEM::Interface & mField
constexpr AssemblyType AT
constexpr IntegrationType IT
#ifndef __CONTACTOPS_HPP__
#define __CONTACTOPS_HPP__
struct CommonData :
public boost::enable_shared_from_this<CommonData> {
static SmartPetscObj<Vec>
constexpr int ghosts[] = {0, 1, 2, 3, 4};
}
const double *t_ptr;
"get array");
t_ptr[4]};
"restore array");
} else {
}
}
return boost::shared_ptr<MatrixDouble>(shared_from_this(),
}
return boost::shared_ptr<MatrixDouble>(shared_from_this(), &
contactDisp);
}
return boost::shared_ptr<MatrixDouble>(shared_from_this(),
}
return boost::shared_ptr<VectorDouble>(shared_from_this(), &
sdfVals);
}
return boost::shared_ptr<MatrixDouble>(shared_from_this(), &
gradsSdf);
}
return boost::shared_ptr<MatrixDouble>(shared_from_this(), &
hessSdf);
}
return boost::shared_ptr<VectorDouble>(shared_from_this(), &
constraintVals);
}
};
#ifdef PYTHON_SDF
struct SDFPython {
SDFPython() = default;
virtual ~SDFPython() = default;
try {
auto main_module = bp::import("__main__");
mainNamespace = main_module.attr("__dict__");
bp::exec_file(py_file.c_str(), mainNamespace, mainNamespace);
sdfFun = mainNamespace["sdf"];
sdfGradFun = mainNamespace["grad_sdf"];
sdfHessFun = mainNamespace["hess_sdf"];
} catch (bp::error_already_set const &) {
PyErr_Print();
}
};
template <typename T>
inline std::vector<T>
py_list_to_std_vector(const boost::python::object &iterable) {
return std::vector<T>(boost::python::stl_input_iterator<T>(iterable),
boost::python::stl_input_iterator<T>());
}
double delta_t,
double t, np::ndarray x, np::ndarray y, np::ndarray z,
np::ndarray tx, np::ndarray ty, np::ndarray tz, int block_id,
) {
try {
sdf = bp::extract<np::ndarray>(
sdfFun(delta_t,
t, x, y, z, tx, ty, tz, block_id));
} catch (bp::error_already_set const &) {
PyErr_Print();
}
}
double delta_t,
double t, np::ndarray x, np::ndarray y, np::ndarray z,
np::ndarray tx, np::ndarray ty, np::ndarray tz, int block_id,
np::ndarray &grad_sdf
) {
try {
sdfGradFun(delta_t,
t, x, y, z, tx, ty, tz, block_id));
} catch (bp::error_already_set const &) {
PyErr_Print();
}
}
double delta_t,
double t, np::ndarray x, np::ndarray y, np::ndarray z,
np::ndarray tx, np::ndarray ty, np::ndarray tz, int block_id,
np::ndarray &hess_sdf
) {
try {
sdfHessFun(delta_t,
t, x, y, z, tx, ty, tz, block_id));
} catch (bp::error_already_set const &) {
PyErr_Print();
}
}
private:
bp::object mainNamespace;
bp::object sdfFun;
bp::object sdfGradFun;
bp::object sdfHessFun;
};
static boost::weak_ptr<SDFPython> sdfPythonWeakPtr;
inline np::ndarray convert_to_numpy(VectorDouble &data, int nb_gauss_pts,
int id) {
auto dtype = np::dtype::get_builtin<double>();
auto size = bp::make_tuple(nb_gauss_pts);
auto stride = bp::make_tuple(3 * sizeof(double));
return (np::from_data(&data[id], dtype, size, stride, bp::object()));
};
#endif
double delta_t,
double t,
int nb_gauss_pts, MatrixDouble &spatial_coords,
MatrixDouble &normals_at_pts, int block_id)>;
double delta_t,
double t,
int nb_gauss_pts, MatrixDouble &spatial_coords,
MatrixDouble &normals_at_pts, int block_id)>;
double delta_t,
double t,
int nb_gauss_pts, MatrixDouble &spatial_coords,
MatrixDouble &normals_at_pts, int block_id)>;
int nb_gauss_pts,
MatrixDouble &m_spatial_coords,
MatrixDouble &m_normals_at_pts,
int block_id) {
#ifdef PYTHON_SDF
if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
VectorDouble v_spatial_coords = m_spatial_coords.data();
VectorDouble v_normal_at_pts = m_normals_at_pts.data();
bp::list python_coords;
bp::list python_normals;
for (int idx = 0; idx < 3; ++idx) {
python_coords.append(
convert_to_numpy(v_spatial_coords, nb_gauss_pts, idx));
python_normals.append(
convert_to_numpy(v_normal_at_pts, nb_gauss_pts, idx));
}
np::ndarray np_sdf = np::empty(bp::make_tuple(nb_gauss_pts),
np::dtype::get_builtin<double>());
bp::extract<np::ndarray>(python_coords[0]),
bp::extract<np::ndarray>(python_coords[1]),
bp::extract<np::ndarray>(python_coords[2]),
bp::extract<np::ndarray>(python_normals[0]),
bp::extract<np::ndarray>(python_normals[1]),
bp::extract<np::ndarray>(python_normals[2]),
block_id, np_sdf),
"Failed python call");
double *sdf_val_ptr = reinterpret_cast<double *>(np_sdf.get_data());
VectorDouble v_sdf;
v_sdf.resize(nb_gauss_pts, false);
for (size_t gg = 0; gg < nb_gauss_pts; ++gg)
v_sdf[gg] = *(sdf_val_ptr + gg);
return v_sdf;
}
#endif
VectorDouble v_sdf;
v_sdf.resize(nb_gauss_pts, false);
auto t_coords = getFTensor1FromPtr<3>(&m_spatial_coords(0, 0));
for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
v_sdf[gg] = -t_coords(2) - 0.1;
++t_coords;
}
return v_sdf;
}
inline MatrixDouble
MatrixDouble &m_spatial_coords,
MatrixDouble &m_normals_at_pts, int block_id) {
#ifdef PYTHON_SDF
if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
VectorDouble v_spatial_coords = m_spatial_coords.data();
VectorDouble v_normal_at_pts = m_normals_at_pts.data();
bp::list python_coords;
bp::list python_normals;
for (int idx = 0; idx < 3; ++idx) {
python_coords.append(
convert_to_numpy(v_spatial_coords, nb_gauss_pts, idx));
python_normals.append(
convert_to_numpy(v_normal_at_pts, nb_gauss_pts, idx));
}
np::ndarray np_grad_sdf = np::empty(bp::make_tuple(nb_gauss_pts, 3),
np::dtype::get_builtin<double>());
delta_t,
t, bp::extract<np::ndarray>(python_coords[0]),
bp::extract<np::ndarray>(python_coords[1]),
bp::extract<np::ndarray>(python_coords[2]),
bp::extract<np::ndarray>(python_normals[0]),
bp::extract<np::ndarray>(python_normals[1]),
bp::extract<np::ndarray>(python_normals[2]), block_id,
np_grad_sdf),
"Failed python call");
double *grad_ptr = reinterpret_cast<double *>(np_grad_sdf.get_data());
MatrixDouble m_grad_sdf;
m_grad_sdf.resize(3, nb_gauss_pts, false);
for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
for (int idx = 0; idx < 3; ++idx)
m_grad_sdf(idx, gg) = *(grad_ptr + (3 * gg + idx));
}
return m_grad_sdf;
}
#endif
MatrixDouble m_grad_sdf;
m_grad_sdf.resize(3, nb_gauss_pts, false);
auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
t_grad_sdf(
i) = t_grad_sdf_set(
i);
++t_grad_sdf;
}
return m_grad_sdf;
}
inline MatrixDouble
MatrixDouble &m_spatial_coords,
MatrixDouble &m_normals_at_pts, int block_id) {
#ifdef PYTHON_SDF
if (auto sdf_ptr = sdfPythonWeakPtr.lock()) {
VectorDouble v_spatial_coords = m_spatial_coords.data();
VectorDouble v_normal_at_pts = m_normals_at_pts.data();
bp::list python_coords;
bp::list python_normals;
for (int idx = 0; idx < 3; ++idx) {
python_coords.append(
convert_to_numpy(v_spatial_coords, nb_gauss_pts, idx));
python_normals.append(
convert_to_numpy(v_normal_at_pts, nb_gauss_pts, idx));
};
np::ndarray np_hess_sdf = np::empty(bp::make_tuple(nb_gauss_pts, 6),
np::dtype::get_builtin<double>());
delta_t,
t, bp::extract<np::ndarray>(python_coords[0]),
bp::extract<np::ndarray>(python_coords[1]),
bp::extract<np::ndarray>(python_coords[2]),
bp::extract<np::ndarray>(python_normals[0]),
bp::extract<np::ndarray>(python_normals[1]),
bp::extract<np::ndarray>(python_normals[2]), block_id,
np_hess_sdf),
"Failed python call");
double *hess_ptr = reinterpret_cast<double *>(np_hess_sdf.get_data());
MatrixDouble m_hess_sdf;
m_hess_sdf.resize(6, nb_gauss_pts, false);
for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
for (int idx = 0; idx < 6; ++idx)
m_hess_sdf(idx, gg) =
*(hess_ptr + (6 * gg + idx));
}
return m_hess_sdf;
}
#endif
MatrixDouble m_hess_sdf;
m_hess_sdf.resize(6, nb_gauss_pts, false);
auto t_hess_sdf = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
for (size_t gg = 0; gg < nb_gauss_pts; ++gg) {
t_hess_sdf(
i,
j) = t_hess_sdf_set(
i,
j);
++t_hess_sdf;
}
return m_hess_sdf;
}
template <int DIM, IntegrationType I, typename BoundaryEleOp>
struct OpAssembleTotalContactTractionImpl;
template <int DIM, IntegrationType I, typename BoundaryEleOp>
struct OpAssembleTotalContactAreaImpl;
template <int DIM, IntegrationType I, typename BoundaryEleOp>
struct OpEvaluateSDFImpl;
template <int DIM, IntegrationType I, typename AssemblyBoundaryEleOp>
struct OpConstrainBoundaryRhsImpl;
template <int DIM, IntegrationType I, typename AssemblyBoundaryEleOp>
struct OpConstrainBoundaryLhs_dUImpl;
template <int DIM, IntegrationType I, typename AssemblyBoundaryEleOp>
struct OpConstrainBoundaryLhs_dTractionImpl;
template <typename T1, typename T2, int DIM1, int DIM2>
size_t nb_gauss_pts) {
MatrixDouble m_spatial_coords(nb_gauss_pts, 3);
m_spatial_coords.clear();
auto t_spatial_coords = getFTensor1FromPtr<3>(&m_spatial_coords(0, 0));
for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
t_spatial_coords(
i) = t_coords(
i) + t_disp(
i);
++t_spatial_coords;
++t_coords;
++t_disp;
}
return m_spatial_coords;
}
template <typename T1, int DIM1>
size_t nb_gauss_pts) {
MatrixDouble m_normals_at_pts(3, nb_gauss_pts);
m_normals_at_pts.clear();
auto t_set_normal = getFTensor1FromMat<3>(m_normals_at_pts);
for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
t_set_normal(
i) = t_normal_at_pts(
i) / t_normal_at_pts.l2();
++t_set_normal;
++t_normal_at_pts;
}
return m_normals_at_pts;
}
template <int DIM, typename BoundaryEleOp>
OpAssembleTotalContactTractionImpl(
boost::shared_ptr<CommonData> common_data_ptr,
double scale = 1,
private:
boost::shared_ptr<CommonData> commonDataPtr;
const double scaleTraction;
bool isAxisymmetric;
};
template <int DIM, typename BoundaryEleOp>
OpAssembleTotalContactAreaImpl(
boost::shared_ptr<CommonData> common_data_ptr,
boost::shared_ptr<Range> contact_range_ptr = nullptr);
private:
boost::shared_ptr<CommonData> commonDataPtr;
bool isAxisymmetric;
boost::shared_ptr<Range> contactRange;
};
template <int DIM, typename BoundaryEleOp>
OpEvaluateSDFImpl(boost::shared_ptr<CommonData> common_data_ptr);
private:
boost::shared_ptr<CommonData> commonDataPtr;
};
template <int DIM, typename AssemblyBoundaryEleOp>
OpConstrainBoundaryRhsImpl(
const std::string
field_name,
boost::shared_ptr<CommonData> common_data_ptr,
private:
boost::shared_ptr<CommonData> commonDataPtr;
bool isAxisymmetric;
};
template <int DIM, typename AssemblyBoundaryEleOp>
OpConstrainBoundaryLhs_dUImpl(const std::string row_field_name,
const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr,
EntitiesFieldData::EntData &col_data);
boost::shared_ptr<CommonData> commonDataPtr;
bool isAxisymmetric;
};
template <int DIM, typename AssemblyBoundaryEleOp>
OpConstrainBoundaryLhs_dTractionImpl(
const std::string row_field_name, const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr,
EntitiesFieldData::EntData &col_data);
private:
boost::shared_ptr<CommonData> commonDataPtr;
bool isAxisymmetric;
};
template <typename BoundaryEleOp> struct ContactIntegrators {
template <int DIM, IntegrationType I>
OpAssembleTotalContactTractionImpl<DIM, I, BoundaryEleOp>;
template <int DIM, IntegrationType I>
OpAssembleTotalContactAreaImpl<DIM, I, BoundaryEleOp>;
template <int DIM, IntegrationType I>
template <AssemblyType A>
struct Assembly {
typename FormsIntegrators<BoundaryEleOp>::template Assembly<A>::OpBase;
template <int DIM, IntegrationType I>
OpConstrainBoundaryRhsImpl<DIM, I, AssemblyBoundaryEleOp>;
template <int DIM, IntegrationType I>
OpConstrainBoundaryLhs_dUImpl<DIM, I, AssemblyBoundaryEleOp>;
template <int DIM, IntegrationType I>
OpConstrainBoundaryLhs_dTractionImpl<DIM, I, AssemblyBoundaryEleOp>;
};
};
inline double sign(
double x) {
constexpr auto eps = std::numeric_limits<float>::epsilon();
return 0;
return 1;
else
return -1;
};
inline double w(
const double sdf,
const double tn) {
}
return (1 - s) / 2;
}
template <int DIM, typename BoundaryEleOp>
OpAssembleTotalContactTractionImpl<DIM, GAUSS, BoundaryEleOp>::
OpAssembleTotalContactTractionImpl(
boost::shared_ptr<CommonData> common_data_ptr,
double scale,
commonDataPtr(common_data_ptr), scaleTraction(
scale),
template <int DIM, typename BoundaryEleOp>
OpAssembleTotalContactTractionImpl<DIM, GAUSS, BoundaryEleOp>::doWork(
int side, EntityType type,
EntData &data) {
auto t_w = BoundaryEleOp::getFTensor0IntegrationWeight();
auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
auto t_coords = BoundaryEleOp::getFTensor1CoordsAtGaussPts();
const auto nb_gauss_pts = BoundaryEleOp::getGaussPts().size2();
for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
double jacobian = 1.;
if (isAxisymmetric) {
jacobian = 2. * M_PI * t_coords(0);
}
const double alpha = t_w * jacobian * BoundaryEleOp::getMeasure();
t_sum_t(
i) += alpha * t_traction(
i);
++t_w;
++t_traction;
++t_coords;
}
t_sum_t(
i) *= scaleTraction;
constexpr int ind[] = {0, 1, 2};
ADD_VALUES);
}
template <int DIM, typename BoundaryEleOp>
OpAssembleTotalContactAreaImpl<DIM, GAUSS, BoundaryEleOp>::
OpAssembleTotalContactAreaImpl(
boost::shared_ptr<Range> contact_range_ptr)
contactRange(contact_range_ptr) {}
template <int DIM, typename BoundaryEleOp>
OpAssembleTotalContactAreaImpl<DIM, GAUSS, BoundaryEleOp>::doWork(
int side, EntityType type,
EntData &data) {
auto fe_type = BoundaryEleOp::getFEType();
const auto fe_ent = BoundaryEleOp::getFEEntityHandle();
if (contactRange->find(fe_ent) != contactRange->end()) {
auto t_w = BoundaryEleOp::getFTensor0IntegrationWeight();
auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
auto t_coords = BoundaryEleOp::getFTensor1CoordsAtGaussPts();
auto t_grad = getFTensor2FromMat<DIM, DIM>(commonDataPtr->contactDispGrad);
auto t_normal_at_pts = BoundaryEleOp::getFTensor1NormalsAtGaussPts();
const auto nb_gauss_pts = BoundaryEleOp::getGaussPts().size2();
BoundaryEleOp::getFTensor1CoordsAtGaussPts(),
getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
BoundaryEleOp::getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
auto t_normal = getFTensor1FromMat<3>(m_normals_at_pts);
auto ts_time = BoundaryEleOp::getTStime();
auto ts_time_step = BoundaryEleOp::getTStimeStep();
int block_id = 0;
auto v_sdf =
surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
m_spatial_coords, m_normals_at_pts, block_id);
auto m_grad_sdf = gradSurfaceDistanceFunction(
ts_time_step, ts_time, nb_gauss_pts, m_spatial_coords, m_normals_at_pts,
block_id);
auto t_sdf = getFTensor0FromVec(v_sdf);
auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
double jacobian = 1.;
if (isAxisymmetric) {
jacobian = 2. * M_PI * t_coords(0);
}
auto tn = -t_traction(
i) * t_grad_sdf(
i);
double alpha = t_w * jacobian;
auto det = determinantTensor(
F);
t_normal_current(
i) = det * (invF(
j,
i) * t_normal_at_pts(
j));
alpha *= sqrt(t_normal_current(
i) * t_normal_current(
i));
if (fe_type == MBTRI) {
alpha /= 2;
}
t_sum_a(0) += alpha;
}
t_sum_a(1) += alpha;
++t_w;
++t_traction;
++t_coords;
++t_sdf;
++t_grad_sdf;
++t_grad;
++t_normal_at_pts;
}
constexpr int ind[] = {3, 4};
ADD_VALUES);
}
}
template <int DIM, typename BoundaryEleOp>
OpEvaluateSDFImpl<DIM, GAUSS, BoundaryEleOp>::OpEvaluateSDFImpl(
boost::shared_ptr<CommonData> common_data_ptr)
commonDataPtr(common_data_ptr) {}
template <int DIM, typename BoundaryEleOp>
OpEvaluateSDFImpl<DIM, GAUSS, BoundaryEleOp>::doWork(int side, EntityType type,
const auto nb_gauss_pts = BoundaryEleOp::getGaussPts().size2();
auto &sdf_vec = commonDataPtr->sdfVals;
auto &grad_mat = commonDataPtr->gradsSdf;
auto &hess_mat = commonDataPtr->hessSdf;
auto &constraint_vec = commonDataPtr->constraintVals;
auto &contactTraction_mat = commonDataPtr->contactTraction;
sdf_vec.resize(nb_gauss_pts, false);
grad_mat.resize(DIM, nb_gauss_pts, false);
hess_mat.resize((DIM * (DIM + 1)) / 2, nb_gauss_pts, false);
constraint_vec.resize(nb_gauss_pts, false);
auto t_traction = getFTensor1FromMat<DIM>(contactTraction_mat);
auto t_sdf = getFTensor0FromVec(sdf_vec);
auto t_grad_sdf = getFTensor1FromMat<DIM>(grad_mat);
auto t_hess_sdf = getFTensor2SymmetricFromMat<DIM>(hess_mat);
auto t_constraint = getFTensor0FromVec(constraint_vec);
auto t_disp = getFTensor1FromMat<DIM>(commonDataPtr->contactDisp);
auto t_coords = BoundaryEleOp::getFTensor1CoordsAtGaussPts();
auto t_normal_at_pts = BoundaryEleOp::getFTensor1NormalsAtGaussPts();
auto ts_time = BoundaryEleOp::getTStime();
auto ts_time_step = BoundaryEleOp::getTStimeStep();
BoundaryEleOp::getFTensor1CoordsAtGaussPts(),
getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
BoundaryEleOp::getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
int block_id = 0;
auto v_sdf =
surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
m_spatial_coords, m_normals_at_pts, block_id);
auto m_grad_sdf =
gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
m_spatial_coords, m_normals_at_pts, block_id);
auto m_hess_sdf =
hessSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
m_spatial_coords, m_normals_at_pts, block_id);
auto t_sdf_v = getFTensor0FromVec(v_sdf);
auto t_grad_sdf_v = getFTensor1FromMat<3>(m_grad_sdf);
auto t_hess_sdf_v = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
auto next = [&]() {
++t_sdf;
++t_sdf_v;
++t_grad_sdf;
++t_grad_sdf_v;
++t_hess_sdf;
++t_hess_sdf_v;
++t_disp;
++t_traction;
++t_constraint;
};
for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
auto tn = -t_traction(
i) * t_grad_sdf_v(
i);
t_sdf = t_sdf_v;
t_grad_sdf(
i) = t_grad_sdf_v(
i);
t_hess_sdf(
i,
j) = t_hess_sdf_v(
i,
j);
next();
}
}
template <int DIM, typename AssemblyBoundaryEleOp>
OpConstrainBoundaryRhsImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::
OpConstrainBoundaryRhsImpl(
const std::string
field_name,
boost::shared_ptr<CommonData> common_data_ptr,
template <int DIM, typename AssemblyBoundaryEleOp>
OpConstrainBoundaryRhsImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::iNtegrate(
EntitiesFieldData::EntData &data) {
const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
auto &nf = AssemblyBoundaryEleOp::locF;
auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
auto t_w = AssemblyBoundaryEleOp::getFTensor0IntegrationWeight();
auto t_disp = getFTensor1FromMat<DIM>(commonDataPtr->contactDisp);
auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
auto t_coords = AssemblyBoundaryEleOp::getFTensor1CoordsAtGaussPts();
size_t nb_base_functions = data.getN().size2() / 3;
auto t_base = data.getFTensor1N<3>();
BoundaryEleOp::getFTensor1CoordsAtGaussPts(),
getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
BoundaryEleOp::getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
auto t_normal = getFTensor1FromMat<3>(m_normals_at_pts);
auto ts_time = AssemblyBoundaryEleOp::getTStime();
auto ts_time_step = AssemblyBoundaryEleOp::getTStimeStep();
int block_id = 0;
auto v_sdf =
surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
m_spatial_coords, m_normals_at_pts, block_id);
auto m_grad_sdf =
gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
m_spatial_coords, m_normals_at_pts, block_id);
auto t_sdf = getFTensor0FromVec(v_sdf);
auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
auto t_nf = getFTensor1FromPtr<DIM>(&nf[0]);
double jacobian = 1.;
if (isAxisymmetric) {
jacobian = 2. * M_PI * t_coords(0);
}
const double alpha = t_w * jacobian * AssemblyBoundaryEleOp::getMeasure();
auto tn = -t_traction(
i) * t_grad_sdf(
i);
t_cP(
i,
j) = (
c * t_grad_sdf(
i)) * t_grad_sdf(
j);
+
c * (t_sdf * t_grad_sdf(
i));
size_t bb = 0;
for (; bb != AssemblyBoundaryEleOp::nbRows / DIM; ++bb) {
const double beta = alpha * (t_base(
i) * t_normal(
i));
t_nf(
i) -= beta * t_rhs(
i);
++t_nf;
++t_base;
}
for (; bb < nb_base_functions; ++bb)
++t_base;
++t_disp;
++t_traction;
++t_coords;
++t_w;
++t_normal;
++t_sdf;
++t_grad_sdf;
}
}
template <int DIM, typename AssemblyBoundaryEleOp>
OpConstrainBoundaryLhs_dUImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::
OpConstrainBoundaryLhs_dUImpl(const std::string row_field_name,
const std::string col_field_name,
boost::shared_ptr<CommonData> common_data_ptr,
AssemblyBoundaryEleOp::sYmm = false;
}
template <int DIM, typename AssemblyBoundaryEleOp>
OpConstrainBoundaryLhs_dUImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::iNtegrate(
EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
auto &locMat = AssemblyBoundaryEleOp::locMat;
auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
auto t_coords = AssemblyBoundaryEleOp::getFTensor1CoordsAtGaussPts();
auto t_w = AssemblyBoundaryEleOp::getFTensor0IntegrationWeight();
auto t_row_base = row_data.getFTensor1N<3>();
size_t nb_face_functions = row_data.getN().size2() / 3;
BoundaryEleOp::getFTensor1CoordsAtGaussPts(),
getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
BoundaryEleOp::getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
auto t_normal = getFTensor1FromMat<3>(m_normals_at_pts);
auto ts_time = AssemblyBoundaryEleOp::getTStime();
auto ts_time_step = AssemblyBoundaryEleOp::getTStimeStep();
int block_id = 0;
auto v_sdf =
surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
m_spatial_coords, m_normals_at_pts, block_id);
auto m_grad_sdf =
gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
m_spatial_coords, m_normals_at_pts, block_id);
auto m_hess_sdf =
hessSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
m_spatial_coords, m_normals_at_pts, block_id);
auto t_sdf = getFTensor0FromVec(v_sdf);
auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
auto t_hess_sdf = getFTensor2SymmetricFromMat<3>(m_hess_sdf);
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
double jacobian = 1.;
if (isAxisymmetric) {
jacobian = 2. * M_PI * t_coords(0);
}
const double alpha = t_w * jacobian * AssemblyBoundaryEleOp::getMeasure();
auto tn = -t_traction(
i) * t_grad_sdf(
i);
t_cP(
i,
j) = (
c * t_grad_sdf(
i)) * t_grad_sdf(
j);
(t_hess_sdf(
i,
j) * (t_grad_sdf(
k) * t_traction(
k)) +
t_grad_sdf(
i) * t_hess_sdf(
k,
j) * t_traction(
k)) +
c * t_sdf * t_hess_sdf(
i,
j);
}
size_t rr = 0;
for (; rr != AssemblyBoundaryEleOp::nbRows / DIM; ++rr) {
auto t_mat = getFTensor2FromArray<DIM, DIM, DIM>(locMat, DIM * rr);
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 / DIM; ++cc) {
const double beta = alpha * row_base * t_col_base;
t_mat(
i,
j) -= beta * t_res_dU(
i,
j);
++t_col_base;
++t_mat;
}
++t_row_base;
}
for (; rr < nb_face_functions; ++rr)
++t_row_base;
++t_traction;
++t_coords;
++t_w;
++t_normal;
++t_sdf;
++t_grad_sdf;
++t_hess_sdf;
}
}
template <int DIM, typename AssemblyBoundaryEleOp>
OpConstrainBoundaryLhs_dTractionImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::
OpConstrainBoundaryLhs_dTractionImpl(
const std::string row_field_name, const std::string col_field_name,
AssemblyBoundaryEleOp::sYmm = false;
}
template <int DIM, typename AssemblyBoundaryEleOp>
OpConstrainBoundaryLhs_dTractionImpl<DIM, GAUSS, AssemblyBoundaryEleOp>::
iNtegrate(EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
const size_t nb_gauss_pts = AssemblyBoundaryEleOp::getGaussPts().size2();
auto &locMat = AssemblyBoundaryEleOp::locMat;
auto t_normal_at_pts = AssemblyBoundaryEleOp::getFTensor1NormalsAtGaussPts();
auto t_traction = getFTensor1FromMat<DIM>(commonDataPtr->contactTraction);
auto t_coords = AssemblyBoundaryEleOp::getFTensor1CoordsAtGaussPts();
auto t_w = AssemblyBoundaryEleOp::getFTensor0IntegrationWeight();
auto t_row_base = row_data.getFTensor1N<3>();
size_t nb_face_functions = row_data.getN().size2() / 3;
BoundaryEleOp::getFTensor1CoordsAtGaussPts(),
getFTensor1FromMat<DIM>(commonDataPtr->contactDisp), nb_gauss_pts);
BoundaryEleOp::getFTensor1NormalsAtGaussPts(), nb_gauss_pts);
auto t_normal = getFTensor1FromMat<3>(m_normals_at_pts);
auto ts_time = AssemblyBoundaryEleOp::getTStime();
auto ts_time_step = AssemblyBoundaryEleOp::getTStimeStep();
int block_id = 0;
auto v_sdf =
surfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
m_spatial_coords, m_normals_at_pts, block_id);
auto m_grad_sdf =
gradSurfaceDistanceFunction(ts_time_step, ts_time, nb_gauss_pts,
m_spatial_coords, m_normals_at_pts, block_id);
auto t_sdf = getFTensor0FromVec(v_sdf);
auto t_grad_sdf = getFTensor1FromMat<3>(m_grad_sdf);
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
double jacobian = 1.;
if (isAxisymmetric) {
jacobian = 2. * M_PI * t_coords(0);
}
const double alpha = t_w * jacobian * AssemblyBoundaryEleOp::getMeasure();
auto tn = -t_traction(
i) * t_grad_sdf(
i);
t_cP(
i,
j) = (
c * t_grad_sdf(
i)) * t_grad_sdf(
j);
size_t rr = 0;
for (; rr != AssemblyBoundaryEleOp::nbRows / DIM; ++rr) {
auto t_mat = getFTensor2FromArray<DIM, DIM, DIM>(locMat, DIM * rr);
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 / 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 * t_res_dt(
i,
j);
++t_col_base;
++t_mat;
}
++t_row_base;
}
for (; rr < nb_face_functions; ++rr)
++t_row_base;
++t_traction;
++t_coords;
++t_w;
++t_normal;
++t_sdf;
++t_grad_sdf;
}
}
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEleOp>
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
using B =
typename FormsIntegrators<DomainEleOp>::template Assembly<
using OpMixDivURhs = typename B::template OpMixDivTimesU<3, DIM, DIM>;
using OpMixDivUCylRhs =
typename B::template OpMixDivTimesU<3, DIM, DIM, CYLINDRICAL>;
using OpMixLambdaGradURhs = typename B::template OpMixTensorTimesGradU<DIM>;
using OpMixUTimesDivLambdaRhs =
typename B::template OpMixVecTimesDivLambda<SPACE_DIM>;
using OpMixUTimesLambdaRhs =
auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
auto mat_grad_ptr = boost::make_shared<MatrixDouble>();
auto div_stress_ptr = boost::make_shared<MatrixDouble>();
auto contact_stress_ptr = boost::make_shared<MatrixDouble>();
else
return 1.;
};
pip.push_back(new OpCalculateVectorFieldValues<DIM>(
u, common_data_ptr->contactDispPtr()));
pip.push_back(
new OpCalculateHVecTensorField<DIM, DIM>(sigma, contact_stress_ptr));
pip.push_back(
new OpCalculateHVecTensorDivergence<DIM, DIM>(sigma, div_stress_ptr));
} else {
pip.push_back(new OpCalculateHVecTensorDivergence<DIM, DIM, CYLINDRICAL>(
sigma, div_stress_ptr));
}
pip.push_back(
new OpMixDivURhs(sigma, common_data_ptr->contactDispPtr(), jacobian));
} else {
pip.push_back(new OpMixDivUCylRhs(sigma, common_data_ptr->contactDispPtr(),
jacobian));
}
pip.push_back(new OpMixLambdaGradURhs(sigma, mat_grad_ptr, jacobian));
pip.push_back(new OpMixUTimesDivLambdaRhs(u, div_stress_ptr, jacobian));
pip.push_back(new OpMixUTimesLambdaRhs(u, contact_stress_ptr, jacobian));
}
template <
typename OpMixLhs>
struct OpMixLhsSide :
public OpMixLhs {
using OpMixLhs::OpMixLhs;
MoFEMErrorCode doWork(
int row_side,
int col_side, EntityType row_type,
EntityType col_type,
EntitiesFieldData::EntData &row_data,
EntitiesFieldData::EntData &col_data) {
auto side_fe_entity = OpMixLhs::getSidePtrFE()->getFEEntityHandle();
auto side_fe_data = OpMixLhs::getSideEntity(row_side, row_type);
if (side_fe_entity == side_fe_data) {
CHKERR OpMixLhs::doWork(row_side, col_side, row_type, col_type, row_data,
col_data);
}
}
};
template <int DIM, AssemblyType A, IntegrationType I, typename DomainEle>
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
std::string fe_domain_name, std::string sigma, std::string u,
std::string geom, ForcesAndSourcesCore::RuleHookFun rule,
using DomainEleOp =
typename DomainEle::UserDataOperator;
auto op_loop_side = new OpLoopSide<DomainEle>(
m_field, fe_domain_name, DIM, Sev::noisy,
boost::make_shared<ForcesAndSourcesCore::UserDataOperator::AdjCache>());
pip.push_back(op_loop_side);
{H1, HDIV}, geom);
using B =
typename FormsIntegrators<DomainEleOp>::template Assembly<
using OpMixDivULhs = typename B::template OpMixDivTimesVec<DIM>;
using OpMixDivUCylLhs =
typename B::template OpMixDivTimesVec<DIM, CYLINDRICAL>;
using OpLambdaGraULhs = typename B::template OpMixTensorTimesGrad<DIM>;
using OpMixDivULhsSide = OpMixLhsSide<OpMixDivULhs>;
using OpMixDivUCylLhsSide = OpMixLhsSide<OpMixDivUCylLhs>;
using OpLambdaGraULhsSide = OpMixLhsSide<OpLambdaGraULhs>;
auto unity = []() { return 1; };
else
return 1.;
};
op_loop_side->getOpPtrVector().push_back(
new OpMixDivULhsSide(sigma, u, unity, jacobian, true));
} else {
op_loop_side->getOpPtrVector().push_back(
new OpMixDivUCylLhsSide(sigma, u, unity, jacobian, true));
}
op_loop_side->getOpPtrVector().push_back(
new OpLambdaGraULhsSide(sigma, u, unity, jacobian, true));
op_loop_side->getSideFEPtr()->getRuleHook = rule;
}
template <int DIM, AssemblyType A, IntegrationType I, typename BoundaryEleOp>
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
using C = ContactIntegrators<BoundaryEleOp>;
auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
pip.push_back(new OpCalculateVectorFieldValues<DIM>(
u, common_data_ptr->contactDispPtr()));
pip.push_back(new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
sigma, common_data_ptr->contactTractionPtr()));
pip.push_back(
new typename C::template Assembly<A>::template OpConstrainBoundaryLhs_dU<
pip.push_back(
new typename C::template Assembly<A>::
template OpConstrainBoundaryLhs_dTraction<DIM, GAUSS>(
}
template <int DIM, AssemblyType A, IntegrationType I, typename BoundaryEleOp>
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
using C = ContactIntegrators<BoundaryEleOp>;
auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
pip.push_back(new OpCalculateVectorFieldValues<DIM>(
u, common_data_ptr->contactDispPtr()));
pip.push_back(new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
sigma, common_data_ptr->contactTractionPtr()));
pip.push_back(
new typename C::template Assembly<A>::template OpConstrainBoundaryRhs<
}
template <int DIM, IntegrationType I, typename BoundaryEleOp>
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
using C = ContactIntegrators<BoundaryEleOp>;
auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
pip.push_back(new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
sigma, common_data_ptr->contactTractionPtr()));
pip.push_back(new typename C::template OpAssembleTotalContactTraction<DIM, I>(
}
template <int DIM, IntegrationType I, typename BoundaryEleOp>
boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pip,
OpLoopSide<SideEle> *op_loop_side, std::string sigma, std::string u,
boost::shared_ptr<Range> contact_range_ptr = nullptr) {
using C = ContactIntegrators<BoundaryEleOp>;
auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
op_loop_side->getOpPtrVector().push_back(
"U", common_data_ptr->contactDispGradPtr()));
if (contact_range_ptr) {
pip.push_back(new OpCalculateVectorFieldValues<DIM>(
u, common_data_ptr->contactDispPtr()));
pip.push_back(new OpCalculateHVecTensorTrace<DIM, BoundaryEleOp>(
sigma, common_data_ptr->contactTractionPtr()));
pip.push_back(op_loop_side);
pip.push_back(new typename C::template OpAssembleTotalContactArea<DIM, I>(
}
}
};
#endif
@ MOFEM_OPERATION_UNSUCCESSFUL
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
Tensor2_Expr< Kronecker_Delta< T >, T, Dim0, Dim1, i, j > kronecker_delta(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
Rank 2.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
MoFEMErrorCode opFactoryDomainRhs(MoFEM::Interface &m_field, std::string block_name, Pip &pip, std::string u, std::string ep, std::string tau)
hess_sdf(t, x, y, z, tx, ty, tz)
grad_sdf(t, x, y, z, tx, ty, tz)
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
FormsIntegrators< BoundaryEleOp >::Assembly< A >::OpBase AssemblyBoundaryEleOp
