15 using SideEle = PipelineManager::ElementsAndOpsByDim<2>::FaceSideEle;
21 using SideEle = PipelineManager::ElementsAndOpsByDim<3>::FaceSideEle;
30 SmartPetscObj<Vec> total_traction);
41 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> ux_scatter,
42 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uy_scatter,
43 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uz_scatter)
48 CHKERR DMoFEMGetInterfacePtr(
dM, &m_field_ptr);
50 createSmartVectorMPI(m_field_ptr->
get_comm(),
53 using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
55 auto get_post_proc_domain_fe = [&]() {
57 boost::make_shared<PostProcEleDomain>(*m_field_ptr,
postProcMesh);
58 auto &pip = post_proc_fe->getOpPtrVector();
60 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
61 auto henky_common_data_ptr = boost::make_shared<HenckyOps::CommonData>();
62 henky_common_data_ptr->matGradPtr = common_data_ptr->mGradPtr();
63 henky_common_data_ptr->matDPtr = common_data_ptr->mDPtr();
65 auto push_domain_ops = [&](
auto &pip) {
68 "Apply base transform");
71 common_data_ptr->mDPtr(), Sev::inform),
74 "U", common_data_ptr->mGradPtr()));
76 new OpCalculateEigenVals<SPACE_DIM>(
"U", henky_common_data_ptr));
78 new OpCalculateLogC<SPACE_DIM>(
"U", henky_common_data_ptr));
80 new OpCalculateLogC_dC<SPACE_DIM>(
"U", henky_common_data_ptr));
82 new OpCalculateHenckyStress<SPACE_DIM>(
"U", henky_common_data_ptr));
84 new OpCalculatePiolaStress<SPACE_DIM>(
"U", henky_common_data_ptr));
85 pip.push_back(
new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(
86 "SIGMA", common_data_ptr->contactStressPtr()));
96 new OpLoopSide<SideEle>(*m_field_ptr,
"dFE",
SPACE_DIM);
98 push_domain_ops(op_loop_side->getOpPtrVector());
99 pip.push_back(op_loop_side);
101 pip.push_back(
new OpCalculateHVecTensorTrace<SPACE_DIM, BoundaryEleOp>(
102 "SIGMA", common_data_ptr->contactTractionPtr()));
104 push_domain_ops(pip);
107 auto u_ptr = boost::make_shared<MatrixDouble>();
108 pip.push_back(
new OpCalculateVectorFieldValues<SPACE_DIM>(
"U", u_ptr));
110 post_proc_fe->getOpPtrVector().push_back(
114 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
122 {
"t", (SPACE_DIM == 3) ? common_data_ptr->contactTractionPtr()
127 {
"SIGMA", common_data_ptr->contactStressPtr()},
129 {
"G", common_data_ptr->mGradPtr()},
131 {
"P2", henky_common_data_ptr->getMatFirstPiolaStress()}
144 auto get_post_proc_bdy_fe = [&]() {
146 boost::make_shared<PostProcEleBdy>(*m_field_ptr,
postProcMesh);
147 auto &pip = post_proc_fe->getOpPtrVector();
149 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
154 pip.push_back(
new OpCalculateVectorFieldValues<SPACE_DIM>(
155 "U", common_data_ptr->contactDispPtr()));
156 pip.push_back(
new OpCalculateHVecTensorTrace<SPACE_DIM, BoundaryEleOp>(
157 "SIGMA", common_data_ptr->contactTractionPtr()));
163 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
167 {{
"U", common_data_ptr->contactDispPtr()},
168 {
"t", common_data_ptr->contactTractionPtr()}},
181 auto get_integrate_traction = [&]() {
182 auto integrate_traction = boost::make_shared<BoundaryEle>(*m_field_ptr);
183 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
185 (AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
186 integrate_traction->getOpPtrVector(), {HDIV} )),
188 integrate_traction->getOpPtrVector().push_back(
189 new OpCalculateHVecTensorTrace<SPACE_DIM, BoundaryEleOp>(
190 "SIGMA", common_data_ptr->contactTractionPtr()));
191 integrate_traction->getOpPtrVector().push_back(
192 new OpAssembleTraction(common_data_ptr, totalTraction));
196 return integrate_traction;
199 postProcDomainFe = get_post_proc_domain_fe();
201 postProcBdyFe = get_post_proc_bdy_fe();
202 integrateTraction = get_integrate_traction();
211 CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
213 auto post_proc = [&]() {
216 auto post_proc_begin =
217 boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(*m_field_ptr,
219 auto post_proc_end = boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
220 *m_field_ptr, postProcMesh);
222 CHKERR DMoFEMPreProcessFiniteElements(dM, post_proc_begin->getFEMethod());
223 if (!postProcBdyFe) {
224 CHKERR DMoFEMLoopFiniteElements(dM,
"bFE", postProcDomainFe);
226 CHKERR DMoFEMLoopFiniteElements(dM,
"dFE", postProcDomainFe);
227 CHKERR DMoFEMLoopFiniteElements(dM,
"bFE", postProcBdyFe);
229 CHKERR DMoFEMPostProcessFiniteElements(dM, post_proc_end->getFEMethod());
231 CHKERR post_proc_end->writeFile(
232 "out_contact_" + boost::lexical_cast<std::string>(sTEP) +
".h5m");
236 auto calculate_traction = [&] {
238 CHKERR VecZeroEntries(totalTraction);
239 CHKERR DMoFEMLoopFiniteElements(dM,
"bFE", integrateTraction);
240 CHKERR VecAssemblyBegin(totalTraction);
241 CHKERR VecAssemblyEnd(totalTraction);
245 auto print_max_min = [&](
auto &tuple,
const std::string msg) {
247 CHKERR VecScatterBegin(std::get<1>(tuple), ts_u, std::get<0>(tuple),
248 INSERT_VALUES, SCATTER_FORWARD);
249 CHKERR VecScatterEnd(std::get<1>(tuple), ts_u, std::get<0>(tuple),
250 INSERT_VALUES, SCATTER_FORWARD);
252 CHKERR VecMax(std::get<0>(tuple), PETSC_NULL, &max);
253 CHKERR VecMin(std::get<0>(tuple), PETSC_NULL, &min);
254 MOFEM_LOG_C(
"CONTACT", Sev::inform,
"%s time %3.4e min %3.4e max %3.4e",
255 msg.c_str(), ts_t, min, max);
259 auto print_traction = [&](
const std::string msg) {
262 CHKERR DMoFEMGetInterfacePtr(dM, &m_field_ptr);
265 CHKERR VecGetArrayRead(totalTraction, &t_ptr);
266 MOFEM_LOG_C(
"CONTACT", Sev::inform,
"%s time %3.4e %3.4e %3.4e %3.4e",
267 msg.c_str(), ts_t, t_ptr[0], t_ptr[1], t_ptr[2]);
268 CHKERR VecRestoreArrayRead(totalTraction, &t_ptr);
274 <<
"Write file at time " << ts_t <<
" write step " << sTEP;
277 CHKERR calculate_traction();
279 CHKERR print_max_min(uXScatter,
"Ux");
280 CHKERR print_max_min(uYScatter,
"Uy");
282 CHKERR print_max_min(uZScatter,
"Uz");
283 CHKERR print_traction(
"Force");
291 SmartPetscObj<DM>
dM;
292 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>>
uXScatter;
293 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>>
uYScatter;
294 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>>
uZScatter;
298 boost::shared_ptr<moab::Core> postProcMesh = boost::make_shared<moab::Core>();
313OpAssembleTraction::OpAssembleTraction(
314 boost::shared_ptr<CommonData> common_data_ptr,
315 SmartPetscObj<Vec> total_traction)
317 commonDataPtr(common_data_ptr), totalTraction(total_traction) {}
327 getFTensor1FromMat<SPACE_DIM>(
commonDataPtr->contactTraction);
330 for (
auto gg = 0; gg != nb_gauss_pts; ++gg) {
332 t_sum_t(
i) += alpha * t_traction(
i);
337 constexpr int ind[] = {0, 1, 2};
#define MOFEM_LOG_C(channel, severity, format,...)
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
@ HDIV
field with continuous normal traction
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MOFEM_LOG(channel, severity)
Log.
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
static constexpr int approx_order
PostProcBrokenMeshInMoabBaseCont< DomainEle > PostProcEleDomain
PipelineManager::ElementsAndOpsByDim< 2 >::FaceSideEle SideEle
PostProcBrokenMeshInMoabBaseCont< BoundaryEle > PostProcEleBdy
PipelineManager::ElementsAndOpsByDim< 3 >::FaceSideEle SideEle
PostProcBrokenMeshInMoabBaseCont< BoundaryEle > PostProcEleBdy
PostProcBrokenMeshInMoabBaseCont< BoundaryEle > PostProcEleDomain
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
default operator for TRI element
double getMeasure()
get measure of element
auto getFTensor0IntegrationWeight()
Get integration weights.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Post post-proc data at points from hash maps.