66 {
68
72
74
76 ->evalFEAtThePoint<DIM>(
81
82 auto process_scalar_field =
83 [](const std::string label,
84 const boost::shared_ptr<VectorDouble> scalarFieldPtr) {
85 if (scalarFieldPtr->size()) {
86 auto t_scalar_holder = getFTensor0FromVec(*scalarFieldPtr);
87
89 << "For " << label << " field, " << t_scalar_holder;
90 }
91 };
92
93 auto process_vector_field =
94 [](const std::string label,
95 const boost::shared_ptr<MatrixDouble> vecFieldPtr) {
96 if (vecFieldPtr->size1()) {
97 auto t_vec_holder = getFTensor1FromMat<DIM>(*vecFieldPtr);
98
99 std::string vec_holder_string = "";
100 for (
int i = 0;
i < DIM;
i++) {
101 vec_holder_string +=
102 " " + boost::lexical_cast<std::string>(t_vec_holder(
i));
103 }
104
106 << "For " << label << " field," << vec_holder_string;
107 }
108 };
109
110 auto process_sym_tensor_field =
111 [](const std::string label,
112 const boost::shared_ptr<MatrixDouble> symTensorFieldPtr) {
113 if (symTensorFieldPtr->size1()) {
114 auto t_sym_tensor_holder =
115 getFTensor2SymmetricFromMat<DIM>(*symTensorFieldPtr);
116
117 std::string sym_tensor_holder_string = "";
118 for (
int i = 0;
i < DIM;
i++) {
119 for (
int j =
i;
j < DIM;
j++) {
120 sym_tensor_holder_string +=
121 ", entry " + boost::lexical_cast<std::string>(
i) +
122 boost::lexical_cast<std::string>(
j) +
" = ";
123 sym_tensor_holder_string +=
124 boost::lexical_cast<std::string>(t_sym_tensor_holder(
i,
j));
125 }
126 }
127
129 << "For " << label << " field" << sym_tensor_holder_string;
130 }
131 };
132
133 auto process_tensor_field =
134 [](const std::string label,
135 const boost::shared_ptr<MatrixDouble> tensorFieldPtr) {
136 if (tensorFieldPtr->size1()) {
137 auto t_tensor_holder =
138 getFTensor2FromMat<DIM, DIM>(*tensorFieldPtr);
139
140 std::string tensor_holder_string = "";
141 for (
int i = 0;
i < DIM;
i++) {
142 for (
int j = 0;
j < DIM;
j++) {
143 tensor_holder_string +=
144 ", entry " + boost::lexical_cast<std::string>(
i) +
145 boost::lexical_cast<std::string>(
j) +
" = ";
146 tensor_holder_string +=
147 boost::lexical_cast<std::string>(t_tensor_holder(
i,
j));
148 }
149 }
150
152 << "For " << label << " field" << tensor_holder_string;
153 }
154 };
155
156 auto process_fields = [](auto &fieldPtr, auto process_field) {
157 for (const auto &pair : fieldPtr) {
158 const std::string &label = pair.first;
159 const boost::shared_ptr<MatrixDouble> &ptr = pair.second;
160
161 if (ptr) {
162 process_field(label, ptr);
163 }
164 }
165 };
166
169 const std::string &label = pair.first;
170 const boost::shared_ptr<VectorDouble> &ptr = pair.second;
171
172 if (ptr) {
173 process_scalar_field(label, ptr);
174 }
175 }
176 }
177
180 const std::string &label = pair.first;
181 const boost::shared_ptr<MatrixDouble> &ptr = pair.second;
182
183 if (ptr) {
184 process_vector_field(label, ptr);
185 }
186 }
187 }
188
191 const std::string &label = pair.first;
192 const boost::shared_ptr<MatrixDouble> &ptr = pair.second;
193
194 if (ptr) {
195 process_sym_tensor_field(label, ptr);
196 }
197 }
198 }
199
202 const std::string &label = pair.first;
203 const boost::shared_ptr<MatrixDouble> &ptr = pair.second;
204
205 if (ptr) {
206 process_tensor_field(label, ptr);
207 }
208 }
209 }
210 }
211
213
214 auto make_vtk = [&]() {
219 "out_plastic_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
220 }
223 getCacheWeakPtr());
225 "out_skin_plastic_" + boost::lexical_cast<std::string>(ts_step) +
226 ".h5m");
227 }
229 };
230
231 auto calculate_reaction = [&]() {
237
238#ifndef NDEBUG
239 auto post_proc_residual = [&](auto dm, auto f_res, auto out_name) {
241 auto post_proc_fe =
242 boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(*m_field_ptr);
243 using OpPPMap = OpPostProcMapInMoab<DIM, DIM>;
244 auto u_vec = boost::make_shared<MatrixDouble>();
245 post_proc_fe->getOpPtrVector().push_back(
246 new OpCalculateVectorFieldValues<DIM>("U", u_vec, f_res));
247 post_proc_fe->getOpPtrVector().push_back(
248
250
251 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
252
253 {},
254
255 {{"RES", u_vec}},
256
257 {}, {})
258
259 );
260
262 post_proc_fe->writeFile("res.h5m");
264 };
265
266 CHKERR post_proc_residual(
dM, r,
"reaction");
267#endif
268
270 };
271
272 auto print_max_min = [&](auto &tuple, const std::string msg) {
274 CHKERR VecScatterBegin(std::get<1>(tuple), ts_u, std::get<0>(tuple),
275 INSERT_VALUES, SCATTER_FORWARD);
276 CHKERR VecScatterEnd(std::get<1>(tuple), ts_u, std::get<0>(tuple),
277 INSERT_VALUES, SCATTER_FORWARD);
278 double max, min;
279 CHKERR VecMax(std::get<0>(tuple), PETSC_NULLPTR, &max);
280 CHKERR VecMin(std::get<0>(tuple), PETSC_NULLPTR, &min);
281 MOFEM_LOG_C(
"PLASTICITY", Sev::inform,
"%s time %3.4e min %3.4e max %3.4e",
282 msg.c_str(), ts_t, min, max);
284 };
285
286 int se = 1;
288 if (!(ts_step % se)) {
290 }
292 CHKERR calculate_reaction();
295 if (DIM == 3)
297
299}
#define MOFEM_LOG_SEVERITY_SYNC(comm, severity)
Synchronise "SYNC" on curtain severity level.
#define MOFEM_LOG_C(channel, severity, format,...)
void simple(double P1[], double P2[], double P3[], double c[], const int N)
#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.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
auto createDMVector(DM dm)
Get smart vector from DM.
PetscErrorCode DMoFEMGetInterfacePtr(DM dm, MoFEM::Interface **m_field_ptr)
Get pointer to MoFEM::Interface.
#define MOFEM_LOG(channel, severity)
Log.
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
PetscBool do_eval_field
Evaluate field.
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
Deprecated interface functions.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.