v0.14.0
Loading...
Searching...
No Matches
ConvectiveMassElement.cpp
Go to the documentation of this file.
1/** \file ConvectiveMassElement.cpp
2 * \brief Operators and data structures for mass and convective mass element
3 * \ingroup convective_mass_elem
4 *
5 */
6
7#include <MoFEM.hpp>
8using namespace MoFEM;
9
11
12#include <adolc/adolc.h>
14#include <DirichletBC.hpp>
17
18#ifndef WITH_ADOL_C
19#error "MoFEM need to be compiled with ADOL-C"
20#endif
21
23 : VolumeElementForcesAndSourcesCore(m_field), A(PETSC_NULL), F(PETSC_NULL) {
24 meshPositionsFieldName = "NoNE";
25
26 auto create_vec = [&]() {
27 constexpr int ghosts[] = {0};
28 if (mField.get_comm_rank() == 0) {
29 return createVectorMPI(mField.get_comm(), 1, 1);
30 } else {
31 return createVectorMPI(mField.get_comm(), 0, 1);
32 }
33 };
34
35 V = create_vec();
36}
37
39
42
43 CHKERR VolumeElementForcesAndSourcesCore::preProcess();
44
45 switch (ts_ctx) {
46 case CTX_TSNONE:
47 CHKERR VecZeroEntries(V);
48 break;
49 default:
50 break;
51 }
52
54}
55
58
59 CHKERR VolumeElementForcesAndSourcesCore::postProcess();
60
61 const double *array;
62 switch (ts_ctx) {
63 case CTX_TSNONE:
64 CHKERR VecAssemblyBegin(V);
65 CHKERR VecAssemblyEnd(V);
66 CHKERR VecSum(V, &eNergy);
67 break;
68 default:
69 break;
70 }
71
73}
74
76 short int tag)
77 : feMassRhs(m_field), feMassLhs(m_field), feMassAuxLhs(m_field),
78 feVelRhs(m_field), feVelLhs(m_field), feTRhs(m_field), feTLhs(m_field),
79 feEnergy(m_field), mField(m_field), tAg(tag) {}
80
82 const std::string field_name,
83 std::vector<VectorDouble> &values_at_gauss_pts,
84 std::vector<MatrixDouble> &gardient_at_gauss_pts)
87 valuesAtGaussPts(values_at_gauss_pts),
88 gradientAtGaussPts(gardient_at_gauss_pts), zeroAtType(MBVERTEX) {}
89
91 int side, EntityType type, EntitiesFieldData::EntData &data) {
93
94 int nb_dofs = data.getFieldData().size();
95 if (nb_dofs == 0) {
97 }
98 int nb_gauss_pts = data.getN().size1();
99 int nb_base_functions = data.getN().size2();
100
101 // initialize
102 // VectorDouble& values = data.getFieldData();
103 valuesAtGaussPts.resize(nb_gauss_pts);
104 gradientAtGaussPts.resize(nb_gauss_pts);
105 for (int gg = 0; gg < nb_gauss_pts; gg++) {
106 valuesAtGaussPts[gg].resize(3);
107 gradientAtGaussPts[gg].resize(3, 3);
108 }
109
110 if (type == zeroAtType) {
111 for (int gg = 0; gg < nb_gauss_pts; gg++) {
112 valuesAtGaussPts[gg].clear();
113 gradientAtGaussPts[gg].clear();
114 }
115 }
116
117 auto base_function = data.getFTensor0N();
118 auto diff_base_functions = data.getFTensor1DiffN<3>();
121
122 for (int gg = 0; gg != nb_gauss_pts; gg++) {
123 auto field_data = data.getFTensor1FieldData<3>();
124 FTensor::Tensor1<double *, 3> values(&valuesAtGaussPts[gg][0],
125 &valuesAtGaussPts[gg][1],
126 &valuesAtGaussPts[gg][2]);
128 &gradientAtGaussPts[gg](0, 0), &gradientAtGaussPts[gg](0, 1),
129 &gradientAtGaussPts[gg](0, 2), &gradientAtGaussPts[gg](1, 0),
130 &gradientAtGaussPts[gg](1, 1), &gradientAtGaussPts[gg](1, 2),
131 &gradientAtGaussPts[gg](2, 0), &gradientAtGaussPts[gg](2, 1),
132 &gradientAtGaussPts[gg](2, 2));
133 int bb = 0;
134 for (; bb != nb_dofs / 3; bb++) {
135 values(i) += base_function * field_data(i);
136 gradient(i, j) += field_data(i) * diff_base_functions(j);
137 ++diff_base_functions;
138 ++base_function;
139 ++field_data;
140 }
141 for (; bb != nb_base_functions; bb++) {
142 ++diff_base_functions;
143 ++base_function;
144 }
145 }
147}
148
150 const std::string field_name, CommonData &common_data)
151 : OpGetDataAtGaussPts(field_name, common_data.dataAtGaussPts[field_name],
152 common_data.gradAtGaussPts[field_name]) {}
153
155 const std::string field_name, BlockData &data, CommonData &common_data,
156 boost::ptr_vector<MethodForForceScaling> &methods_op, int tag,
157 bool jacobian)
160 dAta(data), commonData(common_data), tAg(tag), jAcobian(jacobian),
161 lInear(commonData.lInear), fieldDisp(false), methodsOp(methods_op) {}
162
164 int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data) {
166
167 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
168 dAta.tEts.end()) {
170 }
171
172 // do it only once, no need to repeat this for edges,faces or tets
173 if (row_type != MBVERTEX)
175
176 int nb_dofs = row_data.getIndices().size();
177 if (nb_dofs == 0)
179
180 {
181
182 if (a.size() != 3) {
183 a.resize(3, false);
184 dot_W.resize(3, false);
185 a_res.resize(3, false);
186 g.resize(3, 3, false);
187 G.resize(3, 3, false);
188 h.resize(3, 3, false);
189 H.resize(3, 3, false);
190 invH.resize(3, 3, false);
191 F.resize(3, 3, false);
192 }
193
194 std::fill(dot_W.begin(), dot_W.end(), 0);
195 std::fill(H.data().begin(), H.data().end(), 0);
196 std::fill(invH.data().begin(), invH.data().end(), 0);
197 for (int ii = 0; ii != 3; ii++) {
198 H(ii, ii) = 1;
199 invH(ii, ii) = 1;
200 }
201
202 int nb_gauss_pts = row_data.getN().size1();
203 commonData.valMass.resize(nb_gauss_pts);
204 commonData.jacMassRowPtr.resize(nb_gauss_pts);
205 commonData.jacMass.resize(nb_gauss_pts);
206
207 const std::vector<VectorDouble> &dot_spacial_vel =
209
210 const std::vector<MatrixDouble> &spatial_positions_grad =
212
213 const std::vector<MatrixDouble> &spatial_velocities_grad =
215
216 const std::vector<VectorDouble> &meshpos_vel =
218
219 const std::vector<MatrixDouble> &mesh_positions_gradient =
221
222 int nb_active_vars = 0;
223 for (int gg = 0; gg < nb_gauss_pts; gg++) {
224
225 if (gg == 0) {
226
227 trace_on(tAg);
228
229 for (int nn1 = 0; nn1 < 3; nn1++) { // 0
230 // commonData.dataAtGaussPts["DOT_"+commonData.spatialVelocities]
231 a[nn1] <<= dot_spacial_vel[gg][nn1];
232 nb_active_vars++;
233 }
234 for (int nn1 = 0; nn1 < 3; nn1++) { // 3
235 for (int nn2 = 0; nn2 < 3; nn2++) {
236 // commonData.gradAtGaussPts[commonData.spatialPositions][gg]
237 h(nn1, nn2) <<= spatial_positions_grad[gg](nn1, nn2);
238 if (fieldDisp) {
239 if (nn1 == nn2) {
240 h(nn1, nn2) += 1;
241 }
242 }
243 nb_active_vars++;
244 }
245 }
247 .size() > 0) {
248 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+9=12
249 for (int nn2 = 0; nn2 < 3; nn2++) {
250 // commonData.gradAtGaussPts[commonData.spatialVelocities]
251 g(nn1, nn2) <<= spatial_velocities_grad[gg](nn1, nn2);
252 nb_active_vars++;
253 }
254 }
255 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+9+9=21
256 // commonData.dataAtGaussPts["DOT_"+commonData.meshPositions]
257 dot_W(nn1) <<= meshpos_vel[gg][nn1];
258 nb_active_vars++;
259 }
260 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+9+9+3=24
261 for (int nn2 = 0; nn2 < 3; nn2++) {
262 // commonData.gradAtGaussPts[commonData.meshPositions][gg]
263 H(nn1, nn2) <<= mesh_positions_gradient[gg](nn1, nn2);
264 nb_active_vars++;
265 }
266 }
267 }
268
269 auto a0 = dAta.a0;
271
272 auto t_a_res =
273 FTensor::Tensor1<adouble *, 3>{&a_res[0], &a_res[1], &a_res[2]};
274 auto t_a = FTensor::Tensor1<adouble *, 3>{&a[0], &a[1], &a[2]};
275 auto t_a0 = FTensor::Tensor1<double *, 3>{&a0[0], &a0[1], &a0[2]};
276 auto t_dotW =
277 FTensor::Tensor1<adouble *, 3>{&dot_W[0], &dot_W[1], &dot_W[2]};
280 auto t_invH = getFTensor2FromArray3by3(invH, FTensor::Number<0>(), 0);
283
284 const double rho0 = dAta.rho0;
285
287 CHKERR invertTensor3by3(H, detH, invH);
288
289 t_G(i, j) = t_g(i, k) * t_invH(k, j);
290 t_a_res(i) = t_a(i) - t_a0(i) + t_G(i, j) * t_dotW(j);
291
292 // FIXME: there is error somewhere for nonlinear case
293 // test dam example with -is_linear 0
294 if (!lInear) {
295
296 t_F(i, j) = t_h(i, k) * t_invH(k, j);
297 t_a_res(i) *= rho0 * detH;
298 t_a_res(i) *= determinantTensor3by3(t_F);
299
300 } else {
301
302 t_a_res(i) *= rho0 * detH;
303 }
304
305 // dependant
307 res.resize(3);
308 for (int rr = 0; rr < 3; rr++) {
309 a_res[rr] >>= res[rr];
310 }
311
312 trace_off();
313 }
314
315 active.resize(nb_active_vars);
316 int aa = 0;
317 for (int nn1 = 0; nn1 < 3; nn1++) { // 0
318 active[aa++] = dot_spacial_vel[gg][nn1];
319 }
320 for (int nn1 = 0; nn1 < 3; nn1++) { // 3
321 for (int nn2 = 0; nn2 < 3; nn2++) {
322 if (fieldDisp && nn1 == nn2) {
323 active[aa++] = spatial_positions_grad[gg](nn1, nn2) + 1;
324 } else {
325 active[aa++] = spatial_positions_grad[gg](nn1, nn2);
326 }
327 }
328 }
330 0) {
331 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+9=12
332 for (int nn2 = 0; nn2 < 3; nn2++) {
333 active[aa++] = spatial_velocities_grad[gg](nn1, nn2);
334 }
335 }
336 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+9+9=21
337 active[aa++] = meshpos_vel[gg][nn1];
338 }
339 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+9+9+3=24
340 for (int nn2 = 0; nn2 < 3; nn2++) {
341 active[aa++] = mesh_positions_gradient[gg](nn1, nn2);
342 }
343 }
344 }
345
346 if (!jAcobian) {
348 if (gg > 0) {
349 res.resize(3);
350 int r;
351 r = ::function(tAg, 3, nb_active_vars, &active[0], &res[0]);
352 if (r != 3) { // function is locally analytic
353 SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
354 "ADOL-C function evaluation with error r = %d", r);
355 }
356 }
357 double val = getVolume() * getGaussPts()(3, gg);
358 res *= val;
359 // cout << "my res " << res << endl;
360 } else {
361 commonData.jacMassRowPtr[gg].resize(3);
362 commonData.jacMass[gg].resize(3, nb_active_vars);
363 for (int nn1 = 0; nn1 < 3; nn1++) {
364 (commonData.jacMassRowPtr[gg])[nn1] =
365 &(commonData.jacMass[gg](nn1, 0));
366 }
367 int r;
368 r = jacobian(tAg, 3, nb_active_vars, &active[0],
369 &(commonData.jacMassRowPtr[gg])[0]);
370 if (r != 3) {
371 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
372 "ADOL-C function evaluation with error");
373 }
374 double val = getVolume() * getGaussPts()(3, gg);
375 commonData.jacMass[gg] *= val;
376 }
377 }
378 }
379
381}
382
384 BlockData &data,
385 CommonData &common_data)
388 dAta(data), commonData(common_data) {}
389
392 EntitiesFieldData::EntData &row_data) {
394
395 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
396 dAta.tEts.end()) {
398 }
399 if (row_data.getIndices().size() == 0)
401 int nb_dofs = row_data.getIndices().size();
402
403 auto base = row_data.getFTensor0N();
404 int nb_base_functions = row_data.getN().size2();
405
406 {
407
408 nf.resize(nb_dofs);
409 nf.clear();
410
412
413 for (unsigned int gg = 0; gg < row_data.getN().size1(); gg++) {
414 FTensor::Tensor1<double *, 3> t_nf(&nf[0], &nf[1], &nf[2], 3);
416 &commonData.valMass[gg][1],
417 &commonData.valMass[gg][2]);
418 int dd = 0;
419 for (; dd < nb_dofs / 3; dd++) {
420 t_nf(i) += base * res(i);
421 ++base;
422 ++t_nf;
423 }
424 for (; dd != nb_base_functions; dd++) {
425 ++base;
426 }
427 }
428
429 if ((unsigned int)nb_dofs > 3 * row_data.getN().size2()) {
430 SETERRQ(PETSC_COMM_SELF, 1, "data inconsistency");
431 }
432 CHKERR VecSetValues(getFEMethod()->ts_F, nb_dofs, &row_data.getIndices()[0],
433 &nf[0], ADD_VALUES);
434 }
435
437}
438
440 const std::string vel_field, const std::string field_name, BlockData &data,
441 CommonData &common_data, Range *forcesonlyonentities_ptr)
443 vel_field, field_name,
445 dAta(data), commonData(common_data) {
446 sYmm = false;
447 if (forcesonlyonentities_ptr != NULL) {
448 forcesOnlyOnEntities = *forcesonlyonentities_ptr;
449 }
450}
451
453 EntitiesFieldData::EntData &col_data, int gg) {
455 int nb_col = col_data.getIndices().size();
456 jac.clear();
457 if (!nb_col)
462 FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
463 &jac(1, 0), &jac(1, 1), &jac(1, 2),
464 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
466 &commonData.jacMass[gg](0, 0), &commonData.jacMass[gg](0, 1),
467 &commonData.jacMass[gg](0, 2), &commonData.jacMass[gg](1, 0),
468 &commonData.jacMass[gg](1, 1), &commonData.jacMass[gg](1, 2),
469 &commonData.jacMass[gg](2, 0), &commonData.jacMass[gg](2, 1),
470 &commonData.jacMass[gg](2, 2));
471 double *base_ptr = const_cast<double *>(&col_data.getN(gg)[0]);
472 FTensor::Tensor0<double *> base(base_ptr, 1);
473 if (commonData.dataAtGaussPts["DOT_" + commonData.meshPositions].size() ==
474 0) {
475 for (int dd = 0; dd < nb_col / 3; dd++) {
476 t_jac(i, j) += t_mass1(i, j) * base * getFEMethod()->ts_a;
477 ++base;
478 ++t_jac;
479 }
480 } else {
481 const int s = 3 + 9;
483 // T* d000, T* d001, T* d002,
484 // T* d010, T* d011, T* d012,
485 // T* d020, T* d021, T* d022,
486 // T* d100, T* d101, T* d102,
487 // T* d110, T* d111, T* d112,
488 // T* d120, T* d121, T* d122,
489 // T* d200, T* d201, T* d202,
490 // T* d210, T* d211, T* d212,
491 // T* d220, T* d221, T* d222,
492 &commonData.jacMass[gg](0, s + 0), &commonData.jacMass[gg](0, s + 1),
493 &commonData.jacMass[gg](0, s + 2), &commonData.jacMass[gg](0, s + 3),
494 &commonData.jacMass[gg](0, s + 4), &commonData.jacMass[gg](0, s + 5),
495 &commonData.jacMass[gg](0, s + 6), &commonData.jacMass[gg](0, s + 7),
496 &commonData.jacMass[gg](0, s + 8), &commonData.jacMass[gg](1, s + 0),
497 &commonData.jacMass[gg](1, s + 1), &commonData.jacMass[gg](1, s + 2),
498 &commonData.jacMass[gg](1, s + 3), &commonData.jacMass[gg](1, s + 4),
499 &commonData.jacMass[gg](1, s + 5), &commonData.jacMass[gg](1, s + 6),
500 &commonData.jacMass[gg](1, s + 7), &commonData.jacMass[gg](1, s + 8),
501 &commonData.jacMass[gg](2, s + 0), &commonData.jacMass[gg](2, s + 1),
502 &commonData.jacMass[gg](2, s + 2), &commonData.jacMass[gg](2, s + 3),
503 &commonData.jacMass[gg](2, s + 4), &commonData.jacMass[gg](2, s + 5),
504 &commonData.jacMass[gg](2, s + 6), &commonData.jacMass[gg](2, s + 7),
505 &commonData.jacMass[gg](2, s + 8));
506
507 double *diff_ptr =
508 const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
509 FTensor::Tensor1<double *, 3> diff(diff_ptr, &diff_ptr[1], &diff_ptr[2], 3);
510 for (int dd = 0; dd < nb_col / 3; dd++) {
511 t_jac(i, j) += t_mass1(i, j) * base * getFEMethod()->ts_a;
512 t_jac(i, j) += t_mass3(i, j, k) * diff(k);
513 ++base;
514 ++diff;
515 ++t_jac;
516 }
517 }
519}
520
522 int row_side, int col_side, EntityType row_type, EntityType col_type,
524 EntitiesFieldData::EntData &col_data) {
526
527 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
528 dAta.tEts.end()) {
530 }
531
532 int nb_row = row_data.getIndices().size();
533 int nb_col = col_data.getIndices().size();
534 if (nb_row == 0)
536 if (nb_col == 0)
538
539 auto base = row_data.getFTensor0N();
540 int nb_base_functions = row_data.getN().size2();
541
542 {
543
544 k.resize(nb_row, nb_col);
545 k.clear();
546 jac.resize(3, nb_col);
547
548 for (unsigned int gg = 0; gg < row_data.getN().size1(); gg++) {
549
550 try {
551 CHKERR getJac(col_data, gg);
552 } catch (const std::exception &ex) {
553 std::ostringstream ss;
554 ss << "throw in method: " << ex.what() << std::endl;
555 SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
556 }
557
560
561 {
562 int dd1 = 0;
563 // integrate element stiffness matrix
564 for (; dd1 < nb_row / 3; dd1++) {
566 &jac(0, 0), &jac(0, 1), &jac(0, 2), &jac(1, 0), &jac(1, 1),
567 &jac(1, 2), &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
568 for (int dd2 = 0; dd2 < nb_col / 3; dd2++) {
570 &k(3 * dd1 + 0, 3 * dd2 + 0), &k(3 * dd1 + 0, 3 * dd2 + 1),
571 &k(3 * dd1 + 0, 3 * dd2 + 2), &k(3 * dd1 + 1, 3 * dd2 + 0),
572 &k(3 * dd1 + 1, 3 * dd2 + 1), &k(3 * dd1 + 1, 3 * dd2 + 2),
573 &k(3 * dd1 + 2, 3 * dd2 + 0), &k(3 * dd1 + 2, 3 * dd2 + 1),
574 &k(3 * dd1 + 2, 3 * dd2 + 2));
575 t_k(i, j) += base * t_jac(i, j);
576 ++t_jac;
577 }
578 ++base;
579 // for(int rr1 = 0;rr1<3;rr1++) {
580 // for(int dd2 = 0;dd2<nb_col;dd2++) {
581 // k(3*dd1+rr1,dd2) += row_data.getN()(gg,dd1)*jac(rr1,dd2);
582 // }
583 // }
584 }
585 for (; dd1 != nb_base_functions; dd1++) {
586 ++base;
587 }
588 }
589 }
590
591 if (!forcesOnlyOnEntities.empty()) {
592 VectorInt indices = row_data.getIndices();
593 VectorDofs &dofs = row_data.getFieldDofs();
594 VectorDofs::iterator dit = dofs.begin();
595 for (int ii = 0; dit != dofs.end(); dit++, ii++) {
596 if (forcesOnlyOnEntities.find((*dit)->getEnt()) ==
597 forcesOnlyOnEntities.end()) {
598 indices[ii] = -1;
599 }
600 }
601 CHKERR MatSetValues(getFEMethod()->ts_B, nb_row, &indices[0], nb_col,
602 &col_data.getIndices()[0], &k(0, 0), ADD_VALUES);
603 } else {
604 CHKERR MatSetValues(getFEMethod()->ts_B, nb_row,
605 &row_data.getIndices()[0], nb_col,
606 &col_data.getIndices()[0], &k(0, 0), ADD_VALUES);
607 }
608 }
610}
611
613 const std::string field_name, const std::string col_field, BlockData &data,
614 CommonData &common_data)
615 : OpMassLhs_dM_dv(field_name, col_field, data, common_data) {}
616
618 EntitiesFieldData::EntData &col_data, int gg) {
623 int nb_col = col_data.getIndices().size();
624 jac.clear();
625 FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
626 &jac(1, 0), &jac(1, 1), &jac(1, 2),
627 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
628 const int s = 3;
630 // T* d000, T* d001, T* d002,
631 // T* d010, T* d011, T* d012,
632 // T* d020, T* d021, T* d022,
633 // T* d100, T* d101, T* d102,
634 // T* d110, T* d111, T* d112,
635 // T* d120, T* d121, T* d122,
636 // T* d200, T* d201, T* d202,
637 // T* d210, T* d211, T* d212,
638 // T* d220, T* d221, T* d222,
639 &commonData.jacMass[gg](0, s + 0), &commonData.jacMass[gg](0, s + 1),
640 &commonData.jacMass[gg](0, s + 2), &commonData.jacMass[gg](0, s + 3),
641 &commonData.jacMass[gg](0, s + 4), &commonData.jacMass[gg](0, s + 5),
642 &commonData.jacMass[gg](0, s + 6), &commonData.jacMass[gg](0, s + 7),
643 &commonData.jacMass[gg](0, s + 8), &commonData.jacMass[gg](1, s + 0),
644 &commonData.jacMass[gg](1, s + 1), &commonData.jacMass[gg](1, s + 2),
645 &commonData.jacMass[gg](1, s + 3), &commonData.jacMass[gg](1, s + 4),
646 &commonData.jacMass[gg](1, s + 5), &commonData.jacMass[gg](1, s + 6),
647 &commonData.jacMass[gg](1, s + 7), &commonData.jacMass[gg](1, s + 8),
648 &commonData.jacMass[gg](2, s + 0), &commonData.jacMass[gg](2, s + 1),
649 &commonData.jacMass[gg](2, s + 2), &commonData.jacMass[gg](2, s + 3),
650 &commonData.jacMass[gg](2, s + 4), &commonData.jacMass[gg](2, s + 5),
651 &commonData.jacMass[gg](2, s + 6), &commonData.jacMass[gg](2, s + 7),
652 &commonData.jacMass[gg](2, s + 8));
653 double *diff_ptr =
654 const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
655 FTensor::Tensor1<double *, 3> diff(diff_ptr, &diff_ptr[1], &diff_ptr[2], 3);
656 for (int dd = 0; dd < nb_col / 3; dd++) {
657 t_jac(i, j) += t_mass3(i, j, k) * diff(k);
658 ++diff;
659 ++t_jac;
660 }
662}
663
665 const std::string field_name, const std::string col_field, BlockData &data,
666 CommonData &common_data)
667 : OpMassLhs_dM_dv(field_name, col_field, data, common_data) {}
668
670 EntitiesFieldData::EntData &col_data, int gg) {
672 int nb_col = col_data.getIndices().size();
673 jac.clear();
674 double *base_ptr = const_cast<double *>(&col_data.getN(gg)[0]);
675 FTensor::Tensor0<double *> base(base_ptr, 1);
676 double *diff_ptr =
677 const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
678 FTensor::Tensor1<double *, 3> diff(diff_ptr, &diff_ptr[1], &diff_ptr[2], 3);
679 FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
680 &jac(1, 0), &jac(1, 1), &jac(1, 2),
681 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
682 const int u = 3 + 9 + 9;
684 &commonData.jacMass[gg](0, u + 0), &commonData.jacMass[gg](0, u + 1),
685 &commonData.jacMass[gg](0, u + 2), &commonData.jacMass[gg](1, u + 0),
686 &commonData.jacMass[gg](1, u + 1), &commonData.jacMass[gg](1, u + 2),
687 &commonData.jacMass[gg](2, u + 0), &commonData.jacMass[gg](2, u + 1),
688 &commonData.jacMass[gg](2, u + 2));
689 const int s = 3 + 9 + 9 + 3;
691 // T* d000, T* d001, T* d002,
692 // T* d010, T* d011, T* d012,
693 // T* d020, T* d021, T* d022,
694 // T* d100, T* d101, T* d102,
695 // T* d110, T* d111, T* d112,
696 // T* d120, T* d121, T* d122,
697 // T* d200, T* d201, T* d202,
698 // T* d210, T* d211, T* d212,
699 // T* d220, T* d221, T* d222,
700 &commonData.jacMass[gg](0, s + 0), &commonData.jacMass[gg](0, s + 1),
701 &commonData.jacMass[gg](0, s + 2), &commonData.jacMass[gg](0, s + 3),
702 &commonData.jacMass[gg](0, s + 4), &commonData.jacMass[gg](0, s + 5),
703 &commonData.jacMass[gg](0, s + 6), &commonData.jacMass[gg](0, s + 7),
704 &commonData.jacMass[gg](0, s + 8), &commonData.jacMass[gg](1, s + 0),
705 &commonData.jacMass[gg](1, s + 1), &commonData.jacMass[gg](1, s + 2),
706 &commonData.jacMass[gg](1, s + 3), &commonData.jacMass[gg](1, s + 4),
707 &commonData.jacMass[gg](1, s + 5), &commonData.jacMass[gg](1, s + 6),
708 &commonData.jacMass[gg](1, s + 7), &commonData.jacMass[gg](1, s + 8),
709 &commonData.jacMass[gg](2, s + 0), &commonData.jacMass[gg](2, s + 1),
710 &commonData.jacMass[gg](2, s + 2), &commonData.jacMass[gg](2, s + 3),
711 &commonData.jacMass[gg](2, s + 4), &commonData.jacMass[gg](2, s + 5),
712 &commonData.jacMass[gg](2, s + 6), &commonData.jacMass[gg](2, s + 7),
713 &commonData.jacMass[gg](2, s + 8));
717 for (int dd = 0; dd < nb_col / 3; dd++) {
718 t_jac(i, j) += t_mass1(i, j) * base * getFEMethod()->ts_a;
719 t_jac(i, j) += t_mass3(i, j, k) * diff(k);
720 ++base_ptr;
721 ++diff_ptr;
722 ++t_jac;
723 }
725}
726
728 BlockData &data,
729 CommonData &common_data,
733 dAta(data), commonData(common_data), V(v, true),
734 lInear(commonData.lInear) {}
735
738 EntitiesFieldData::EntData &row_data) {
740
741 if (row_type != MBVERTEX) {
743 }
744 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
745 dAta.tEts.end()) {
747 }
748
749 {
750 double energy = 0;
751 for (unsigned int gg = 0; gg < row_data.getN().size1(); gg++) {
752 double val = getVolume() * getGaussPts()(3, gg);
753 double rho0 = dAta.rho0;
754 double rho;
755 if (lInear) {
756 rho = rho0;
757 } else {
758 h.resize(3, 3);
759 noalias(h) =
762 .size() > 0) {
763 H.resize(3, 3);
764 noalias(H) =
766 auto detH = determinantTensor3by3(H);
767 invH.resize(3, 3);
768 CHKERR invertTensor3by3(H, detH, invH);
769 F.resize(3, 3);
770 noalias(F) = prod(h, invH);
771 } else {
772 F.resize(3, 3);
773 noalias(F) = h;
774 }
775 double detF = determinantTensor3by3(F);
776 rho = detF * rho0;
777 }
778 v.resize(3);
780 energy += 0.5 * (rho * val) * inner_prod(v, v);
781 }
782 CHKERR VecSetValue(V, 0, energy, ADD_VALUES);
783 }
784
786}
787
789 const std::string field_name, BlockData &data, CommonData &common_data,
790 int tag, bool jacobian)
793 dAta(data), commonData(common_data), tAg(tag), jAcobian(jacobian),
794 fieldDisp(false) {}
795
797 int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data) {
799
800 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
801 dAta.tEts.end()) {
803 }
804
805 // do it only once, no need to repeat this for edges,faces or tets
806 if (row_type != MBVERTEX)
808
809 int nb_dofs = row_data.getIndices().size();
810 if (nb_dofs == 0)
812
813 {
814
815 v.resize(3);
816 dot_w.resize(3);
817 h.resize(3, 3);
818 h.clear();
819 F.resize(3, 3);
820 dot_W.resize(3);
821 dot_W.clear();
822 H.resize(3, 3);
823 H.clear();
824 invH.resize(3, 3);
825 invH.clear();
826 dot_u.resize(3);
827 for (int dd = 0; dd < 3; dd++) {
828 H(dd, dd) = 1;
829 invH(dd, dd) = 1;
830 }
831
832 a_res.resize(3);
833 int nb_gauss_pts = row_data.getN().size1();
834 commonData.valVel.resize(nb_gauss_pts);
835 commonData.jacVelRowPtr.resize(nb_gauss_pts);
836 commonData.jacVel.resize(nb_gauss_pts);
837
838 int nb_active_vars = 0;
839 for (int gg = 0; gg < nb_gauss_pts; gg++) {
840
841 if (gg == 0) {
842
843 trace_on(tAg);
844
845 for (int nn1 = 0; nn1 < 3; nn1++) { // 0
846 v[nn1] <<=
848 nb_active_vars++;
849 }
850 for (int nn1 = 0; nn1 < 3; nn1++) { // 3
851 dot_w[nn1] <<=
853 [gg][nn1];
854 nb_active_vars++;
855 }
857 .size() > 0) {
858 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3 = 6
859 for (int nn2 = 0; nn2 < 3; nn2++) {
860 h(nn1, nn2) <<=
862 nn1, nn2);
863 if (fieldDisp) {
864 if (nn1 == nn2) {
865 h(nn1, nn2) += 1;
866 }
867 }
868 nb_active_vars++;
869 }
870 }
871 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9
872 dot_W[nn1] <<=
874 .dataAtGaussPts["DOT_" + commonData.meshPositions][gg][nn1];
875 nb_active_vars++;
876 }
877 }
879 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9+3
880 for (int nn2 = 0; nn2 < 3; nn2++) {
881 H(nn1, nn2) <<=
883 nn2);
884 nb_active_vars++;
885 }
886 }
887 }
888 detH = 1;
889
893
897 auto t_invH = getFTensor2FromArray3by3(invH, FTensor::Number<0>(), 0);
898 auto t_dot_u =
899 FTensor::Tensor1<adouble *, 3>{&dot_u[0], &dot_u[1], &dot_u[2]};
900 auto t_dot_w =
901 FTensor::Tensor1<adouble *, 3>{&dot_w[0], &dot_w[1], &dot_w[2]};
902 auto t_dot_W =
903 FTensor::Tensor1<adouble *, 3>{&dot_W[0], &dot_W[1], &dot_W[2]};
904 auto t_v = FTensor::Tensor1<adouble *, 3>{&v[0], &v[1], &v[2]};
905 auto t_a_res =
906 FTensor::Tensor1<adouble *, 3>{&a_res[0], &a_res[1], &a_res[2]};
907
909 detH = determinantTensor3by3(H);
910 CHKERR invertTensor3by3(H, detH, invH);
911 t_F(i, j) = t_h(i, k) * t_invH(k, j);
912 } else {
913 t_F(i, j) = t_h(i, j);
914 }
915
916 t_dot_u(i) = t_dot_w(i) + t_F(i, j) * t_dot_W(j);
917 t_a_res(i) = t_v(i) - t_dot_u(i);
918 t_a_res(i) *= detH;
919
920 // dependant
921 VectorDouble &res = commonData.valVel[gg];
922 res.resize(3);
923 for (int rr = 0; rr < 3; rr++) {
924 a_res[rr] >>= res[rr];
925 }
926 trace_off();
927 }
928
929 active.resize(nb_active_vars);
930 int aa = 0;
931 for (int nn1 = 0; nn1 < 3; nn1++) {
932 active[aa++] =
934 }
935 for (int nn1 = 0; nn1 < 3; nn1++) {
936 active[aa++] =
938 .dataAtGaussPts["DOT_" + commonData.spatialPositions][gg][nn1];
939 }
941 0) {
942 for (int nn1 = 0; nn1 < 3; nn1++) {
943 for (int nn2 = 0; nn2 < 3; nn2++) {
944 if (fieldDisp && nn1 == nn2) {
945 active[aa++] =
947 nn1, nn2) +
948 1;
949 } else {
950 active[aa++] =
952 nn1, nn2);
953 }
954 }
955 }
956 for (int nn1 = 0; nn1 < 3; nn1++) {
957 active[aa++] =
959 .dataAtGaussPts["DOT_" + commonData.meshPositions][gg][nn1];
960 }
961 }
963 for (int nn1 = 0; nn1 < 3; nn1++) {
964 for (int nn2 = 0; nn2 < 3; nn2++) {
965 active[aa++] =
967 nn2);
968 }
969 }
970 }
971
972 if (!jAcobian) {
973 VectorDouble &res = commonData.valVel[gg];
974 if (gg > 0) {
975 res.resize(3);
976 int r;
977 r = ::function(tAg, 3, nb_active_vars, &active[0], &res[0]);
978 if (r != 3) {
979 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
980 "ADOL-C function evaluation with error");
981 }
982 }
983 double val = getVolume() * getGaussPts()(3, gg);
984 res *= val;
985 } else {
986 commonData.jacVelRowPtr[gg].resize(3);
987 commonData.jacVel[gg].resize(3, nb_active_vars);
988 for (int nn1 = 0; nn1 < 3; nn1++) {
989 (commonData.jacVelRowPtr[gg])[nn1] = &(commonData.jacVel[gg](nn1, 0));
990 }
991 int r;
992 r = jacobian(tAg, 3, nb_active_vars, &active[0],
993 &(commonData.jacVelRowPtr[gg])[0]);
994 if (r != 3) {
995 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
996 "ADOL-C function evaluation with error");
997 }
998 double val = getVolume() * getGaussPts()(3, gg);
999 commonData.jacVel[gg] *= val;
1000 // std::cerr << gg << " : " << commonData.jacVel[gg] << std::endl;
1001 }
1002 }
1003 }
1005}
1006
1008 const std::string field_name, BlockData &data, CommonData &common_data)
1011 dAta(data), commonData(common_data) {}
1012
1014 int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data) {
1016
1017 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
1018 dAta.tEts.end()) {
1020 }
1021 int nb_dofs = row_data.getIndices().size();
1022 if (nb_dofs == 0)
1024
1025 auto base = row_data.getFTensor0N();
1026 int nb_base_functions = row_data.getN().size2();
1028
1029 {
1030
1031 nf.resize(nb_dofs);
1032 nf.clear();
1033
1034 for (unsigned int gg = 0; gg < row_data.getN().size1(); gg++) {
1035 FTensor::Tensor1<double *, 3> t_nf(&nf[0], &nf[1], &nf[2], 3);
1037 &commonData.valVel[gg][1],
1038 &commonData.valVel[gg][2]);
1039 int dd = 0;
1040 for (; dd < nb_dofs / 3; dd++) {
1041 t_nf(i) += base * res(i);
1042 ++base;
1043 ++t_nf;
1044 }
1045 for (; dd != nb_base_functions; dd++) {
1046 ++base;
1047 }
1048 }
1049
1050 if (row_data.getIndices().size() > 3 * row_data.getN().size2()) {
1051 SETERRQ(PETSC_COMM_SELF, 1, "data inconsistency");
1052 }
1053 CHKERR VecSetValues(getFEMethod()->ts_F, row_data.getIndices().size(),
1054 &row_data.getIndices()[0], &nf[0], ADD_VALUES);
1055 }
1057}
1058
1060 const std::string vel_field, const std::string field_name, BlockData &data,
1061 CommonData &common_data)
1062 : OpMassLhs_dM_dv(vel_field, field_name, data, common_data) {}
1063
1065 EntitiesFieldData::EntData &col_data, int gg) {
1067 int nb_col = col_data.getIndices().size();
1068 jac.clear();
1069 if (!nb_col)
1071 double *base_ptr = const_cast<double *>(&col_data.getN(gg)[0]);
1072 FTensor::Tensor0<double *> base(base_ptr, 1);
1073 FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
1074 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1075 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1077 &commonData.jacVel[gg](0, 0), &commonData.jacVel[gg](0, 1),
1078 &commonData.jacVel[gg](0, 2), &commonData.jacVel[gg](1, 0),
1079 &commonData.jacVel[gg](1, 1), &commonData.jacVel[gg](1, 2),
1080 &commonData.jacVel[gg](2, 0), &commonData.jacVel[gg](2, 1),
1081 &commonData.jacVel[gg](2, 2));
1084 for (int dd = 0; dd < nb_col / 3; dd++) {
1085 t_jac(i, j) += t_mass1(i, j) * base;
1086 ++base_ptr;
1087 ++t_jac;
1088 }
1089
1091}
1092
1094 const std::string vel_field, const std::string field_name, BlockData &data,
1095 CommonData &common_data)
1096 : OpVelocityLhs_dV_dv(vel_field, field_name, data, common_data) {}
1097
1099 EntitiesFieldData::EntData &col_data, int gg) {
1101 int nb_col = col_data.getIndices().size();
1102 jac.clear();
1103 if (!nb_col)
1105 double *base_ptr = const_cast<double *>(&col_data.getN(gg)[0]);
1106 FTensor::Tensor0<double *> base(base_ptr, 1);
1107 FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
1108 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1109 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1110 const int u = 3;
1112 &commonData.jacVel[gg](0, u + 0), &commonData.jacVel[gg](0, u + 1),
1113 &commonData.jacVel[gg](0, u + 2), &commonData.jacVel[gg](1, u + 0),
1114 &commonData.jacVel[gg](1, u + 1), &commonData.jacVel[gg](1, u + 2),
1115 &commonData.jacVel[gg](2, u + 0), &commonData.jacVel[gg](2, u + 1),
1116 &commonData.jacVel[gg](2, u + 2));
1120 if (commonData.dataAtGaussPts["DOT_" + commonData.meshPositions].size() ==
1121 0) {
1122
1123 for (int dd = 0; dd < nb_col / 3; dd++) {
1124 t_jac(i, j) += t_mass1(i, j) * base * getFEMethod()->ts_a;
1125 ++base_ptr;
1126 ++t_jac;
1127 }
1128 } else {
1129 double *diff_ptr =
1130 const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
1131 FTensor::Tensor1<double *, 3> diff(diff_ptr, &diff_ptr[1], &diff_ptr[2], 3);
1132 const int s = 3 + 3;
1134 // T* d000, T* d001, T* d002,
1135 // T* d010, T* d011, T* d012,
1136 // T* d020, T* d021, T* d022,
1137 // T* d100, T* d101, T* d102,
1138 // T* d110, T* d111, T* d112,
1139 // T* d120, T* d121, T* d122,
1140 // T* d200, T* d201, T* d202,
1141 // T* d210, T* d211, T* d212,
1142 // T* d220, T* d221, T* d222,
1143 &commonData.jacVel[gg](0, s + 0), &commonData.jacVel[gg](0, s + 1),
1144 &commonData.jacVel[gg](0, s + 2), &commonData.jacVel[gg](0, s + 3),
1145 &commonData.jacVel[gg](0, s + 4), &commonData.jacVel[gg](0, s + 5),
1146 &commonData.jacVel[gg](0, s + 6), &commonData.jacVel[gg](0, s + 7),
1147 &commonData.jacVel[gg](0, s + 8), &commonData.jacVel[gg](1, s + 0),
1148 &commonData.jacVel[gg](1, s + 1), &commonData.jacVel[gg](1, s + 2),
1149 &commonData.jacVel[gg](1, s + 3), &commonData.jacVel[gg](1, s + 4),
1150 &commonData.jacVel[gg](1, s + 5), &commonData.jacVel[gg](1, s + 6),
1151 &commonData.jacVel[gg](1, s + 7), &commonData.jacVel[gg](1, s + 8),
1152 &commonData.jacVel[gg](2, s + 0), &commonData.jacVel[gg](2, s + 1),
1153 &commonData.jacVel[gg](2, s + 2), &commonData.jacVel[gg](2, s + 3),
1154 &commonData.jacVel[gg](2, s + 4), &commonData.jacVel[gg](2, s + 5),
1155 &commonData.jacVel[gg](2, s + 6), &commonData.jacVel[gg](2, s + 7),
1156 &commonData.jacVel[gg](2, s + 8));
1157 for (int dd = 0; dd < nb_col / 3; dd++) {
1158 t_jac(i, j) += t_mass1(i, j) * base * getFEMethod()->ts_a;
1159 t_jac(i, j) += t_mass3(i, j, k) * diff(k);
1160 ++base_ptr;
1161 ++diff_ptr;
1162 ++t_jac;
1163 }
1164 }
1166}
1167
1169 const std::string vel_field, const std::string field_name, BlockData &data,
1170 CommonData &common_data)
1171 : OpVelocityLhs_dV_dv(vel_field, field_name, data, common_data) {}
1172
1174 EntitiesFieldData::EntData &col_data, int gg) {
1176 int nb_col = col_data.getIndices().size();
1177 jac.clear();
1178 if (!nb_col)
1180 double *base_ptr = const_cast<double *>(&col_data.getN(gg)[0]);
1181 FTensor::Tensor0<double *> base(base_ptr, 1);
1182 double *diff_ptr =
1183 const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
1184 FTensor::Tensor1<double *, 3> diff(diff_ptr, &diff_ptr[1], &diff_ptr[2], 3);
1185 FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
1186 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1187 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1188 const int u = 3 + 3 + 9;
1190 &commonData.jacVel[gg](0, u + 0), &commonData.jacVel[gg](0, u + 1),
1191 &commonData.jacVel[gg](0, u + 2), &commonData.jacVel[gg](1, u + 0),
1192 &commonData.jacVel[gg](1, u + 1), &commonData.jacVel[gg](1, u + 2),
1193 &commonData.jacVel[gg](2, u + 0), &commonData.jacVel[gg](2, u + 1),
1194 &commonData.jacVel[gg](2, u + 2));
1195 const int s = 3 + 3 + 9 + 3;
1197 // T* d000, T* d001, T* d002,
1198 // T* d010, T* d011, T* d012,
1199 // T* d020, T* d021, T* d022,
1200 // T* d100, T* d101, T* d102,
1201 // T* d110, T* d111, T* d112,
1202 // T* d120, T* d121, T* d122,
1203 // T* d200, T* d201, T* d202,
1204 // T* d210, T* d211, T* d212,
1205 // T* d220, T* d221, T* d222,
1206 &commonData.jacVel[gg](0, s + 0), &commonData.jacVel[gg](0, s + 1),
1207 &commonData.jacVel[gg](0, s + 2), &commonData.jacVel[gg](0, s + 3),
1208 &commonData.jacVel[gg](0, s + 4), &commonData.jacVel[gg](0, s + 5),
1209 &commonData.jacVel[gg](0, s + 6), &commonData.jacVel[gg](0, s + 7),
1210 &commonData.jacVel[gg](0, s + 8), &commonData.jacVel[gg](1, s + 0),
1211 &commonData.jacVel[gg](1, s + 1), &commonData.jacVel[gg](1, s + 2),
1212 &commonData.jacVel[gg](1, s + 3), &commonData.jacVel[gg](1, s + 4),
1213 &commonData.jacVel[gg](1, s + 5), &commonData.jacVel[gg](1, s + 6),
1214 &commonData.jacVel[gg](1, s + 7), &commonData.jacVel[gg](1, s + 8),
1215 &commonData.jacVel[gg](2, s + 0), &commonData.jacVel[gg](2, s + 1),
1216 &commonData.jacVel[gg](2, s + 2), &commonData.jacVel[gg](2, s + 3),
1217 &commonData.jacVel[gg](2, s + 4), &commonData.jacVel[gg](2, s + 5),
1218 &commonData.jacVel[gg](2, s + 6), &commonData.jacVel[gg](2, s + 7),
1219 &commonData.jacVel[gg](2, s + 8));
1223 for (int dd = 0; dd < nb_col / 3; dd++) {
1224 t_jac(i, j) += t_mass1(i, j) * base * getFEMethod()->ts_a;
1225 t_jac(i, j) += t_mass3(i, j, k) * diff(k);
1226 ++base_ptr;
1227 ++diff_ptr;
1228 ++t_jac;
1229 }
1231}
1232
1235 BlockData &data,
1236 CommonData &common_data, int tag,
1237 bool jacobian)
1240 dAta(data), commonData(common_data), tAg(tag), jAcobian(jacobian),
1241 fieldDisp(false) {}
1242
1245 int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data) {
1247
1248 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
1249 dAta.tEts.end()) {
1251 }
1252
1253 // do it only once, no need to repeat this for edges,faces or tets
1254 if (row_type != MBVERTEX)
1256
1257 int nb_dofs = row_data.getIndices().size();
1258 if (nb_dofs == 0)
1260
1261 try {
1262
1263 a.resize(3);
1264 v.resize(3);
1265 g.resize(3, 3);
1266 G.resize(3, 3);
1267 h.resize(3, 3);
1268 F.resize(3, 3);
1269 H.resize(3, 3);
1270 H.clear();
1271 invH.resize(3, 3);
1272 invH.clear();
1273 for (int dd = 0; dd < 3; dd++) {
1274 H(dd, dd) = 1;
1275 invH(dd, dd) = 1;
1276 }
1277
1278 int nb_gauss_pts = row_data.getN().size1();
1279 commonData.valT.resize(nb_gauss_pts);
1280 commonData.jacTRowPtr.resize(nb_gauss_pts);
1281 commonData.jacT.resize(nb_gauss_pts);
1282
1283 int nb_active_vars = 0;
1284 for (int gg = 0; gg < nb_gauss_pts; gg++) {
1285
1286 if (gg == 0) {
1287
1288 trace_on(tAg);
1289
1290 for (int nn1 = 0; nn1 < 3; nn1++) { // 0
1291 a[nn1] <<=
1293 [gg][nn1];
1294 nb_active_vars++;
1295 }
1296
1297 for (int nn1 = 0; nn1 < 3; nn1++) { // 3
1298 v[nn1] <<=
1300 nb_active_vars++;
1301 }
1302 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3
1303 for (int nn2 = 0; nn2 < 3; nn2++) {
1304 g(nn1, nn2) <<=
1306 nn1, nn2);
1307 nb_active_vars++;
1308 }
1309 }
1310 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9
1311 for (int nn2 = 0; nn2 < 3; nn2++) {
1312 h(nn1, nn2) <<=
1314 nn2);
1315 nb_active_vars++;
1316 if (fieldDisp) {
1317 if (nn1 == nn2) {
1318 h(nn1, nn2) += 1;
1319 }
1320 }
1321 }
1322 }
1324 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9+9
1325 for (int nn2 = 0; nn2 < 3; nn2++) {
1326 H(nn1, nn2) <<=
1328 nn2);
1329 nb_active_vars++;
1330 }
1331 }
1332 }
1333 adouble detH;
1334 detH = 1;
1336 detH = determinantTensor3by3(H);
1337 }
1338 CHKERR invertTensor3by3(H, detH, invH);
1339
1343
1344 a_T.resize(3);
1345
1347 auto t_invH = getFTensor2FromArray3by3(invH, FTensor::Number<0>(), 0);
1350 auto t_G = getFTensor2FromArray3by3(G, FTensor::Number<0>(), 0);
1351
1352 auto t_a = FTensor::Tensor1<adouble *, 3>{&a[0], &a[1], &a[2]};
1353 auto t_v = FTensor::Tensor1<adouble *, 3>{&v[0], &v[1], &v[2]};
1354 auto t_a_T = FTensor::Tensor1<adouble *, 3>{&a_T[0], &a_T[1], &a_T[2]};
1355
1356 t_F(i, j) = t_h(i, k) * t_invH(k, j);
1357 t_G(i, j) = t_g(i, k) * t_invH(k, j);
1358 t_a_T(i) = t_F(k, i) * t_a(k) + t_G(k, i) * t_v(k);
1359 const auto rho0 = dAta.rho0;
1360 t_a_T(i) = -rho0 * detH;
1361
1362 commonData.valT[gg].resize(3);
1363 for (int nn = 0; nn < 3; nn++) {
1364 a_T[nn] >>= (commonData.valT[gg])[nn];
1365 }
1366 trace_off();
1367 }
1368
1369 active.resize(nb_active_vars);
1370 int aa = 0;
1371 for (int nn1 = 0; nn1 < 3; nn1++) { // 0
1372 active[aa++] =
1374 .dataAtGaussPts["DOT_" + commonData.spatialVelocities][gg][nn1];
1375 }
1376
1377 for (int nn1 = 0; nn1 < 3; nn1++) { // 3
1378 active[aa++] =
1380 }
1381 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3
1382 for (int nn2 = 0; nn2 < 3; nn2++) {
1383 active[aa++] =
1385 nn2);
1386 }
1387 }
1388 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9
1389 for (int nn2 = 0; nn2 < 3; nn2++) {
1390 if (fieldDisp && nn1 == nn2) {
1391 active[aa++] =
1393 nn1, nn2) +
1394 1;
1395 } else {
1396 active[aa++] =
1398 nn2);
1399 }
1400 }
1401 }
1403 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9+9
1404 for (int nn2 = 0; nn2 < 3; nn2++) {
1405 active[aa++] =
1407 nn2);
1408 }
1409 }
1410 }
1411
1412 if (!jAcobian) {
1413 VectorDouble &res = commonData.valT[gg];
1414 if (gg > 0) {
1415 res.resize(3);
1416 int r;
1417 r = ::function(tAg, 3, nb_active_vars, &active[0], &res[0]);
1418 if (r != 3) { // function is locally analytic
1419 SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1420 "ADOL-C function evaluation with error r = %d", r);
1421 }
1422 }
1423 double val = getVolume() * getGaussPts()(3, gg);
1424 res *= val;
1425 } else {
1426 commonData.jacTRowPtr[gg].resize(3);
1427 commonData.jacT[gg].resize(3, nb_active_vars);
1428 for (int nn1 = 0; nn1 < 3; nn1++) {
1429 (commonData.jacTRowPtr[gg])[nn1] = &(commonData.jacT[gg](nn1, 0));
1430 }
1431 int r;
1432 r = jacobian(tAg, 3, nb_active_vars, &active[0],
1433 &(commonData.jacTRowPtr[gg])[0]);
1434 if (r != 3) {
1435 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1436 "ADOL-C function evaluation with error");
1437 }
1438 double val = getVolume() * getGaussPts()(3, gg);
1439 commonData.jacT[gg] *= val;
1440 }
1441 }
1442
1443 } catch (const std::exception &ex) {
1444 std::ostringstream ss;
1445 ss << "throw in method: " << ex.what() << std::endl;
1446 SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
1447 }
1448
1450}
1451
1454 BlockData &data,
1455 CommonData &common_data,
1456 Range *forcesonlyonentities_ptr)
1459 dAta(data), commonData(common_data) {
1460 if (forcesonlyonentities_ptr != NULL) {
1461 forcesOnlyOnEntities = *forcesonlyonentities_ptr;
1462 }
1463}
1464
1467 int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data) {
1469
1470 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
1471 dAta.tEts.end()) {
1473 }
1474 int nb_dofs = row_data.getIndices().size();
1475 if (nb_dofs == 0)
1477
1478 try {
1479
1480 nf.resize(nb_dofs);
1481 nf.clear();
1482
1483 auto base = row_data.getFTensor0N();
1484 int nb_base_functions = row_data.getN().size2();
1486
1487 for (unsigned int gg = 0; gg < row_data.getN().size1(); gg++) {
1488 FTensor::Tensor1<double *, 3> t_nf(&nf[0], &nf[1], &nf[2], 3);
1490 &commonData.valT[gg][1],
1491 &commonData.valT[gg][2]);
1492 int dd = 0;
1493 for (; dd < nb_dofs / 3; dd++) {
1494 t_nf(i) += base * res(i);
1495 ++base;
1496 ++t_nf;
1497 }
1498 for (; dd != nb_base_functions; dd++) {
1499 ++base;
1500 }
1501 }
1502
1503 if (row_data.getIndices().size() > 3 * row_data.getN().size2()) {
1504 SETERRQ(PETSC_COMM_SELF, 1, "data inconsistency");
1505 }
1506 if (!forcesOnlyOnEntities.empty()) {
1507 VectorInt indices = row_data.getIndices();
1508 VectorDofs &dofs = row_data.getFieldDofs();
1509 VectorDofs::iterator dit = dofs.begin();
1510 for (int ii = 0; dit != dofs.end(); dit++, ii++) {
1511 if (forcesOnlyOnEntities.find((*dit)->getEnt()) ==
1512 forcesOnlyOnEntities.end()) {
1513 // std::cerr << **dit << std::endl;
1514 indices[ii] = -1;
1515 }
1516 }
1517 // std::cerr << indices << std::endl;
1518 CHKERR VecSetValues(getFEMethod()->ts_F, indices.size(), &indices[0],
1519 &nf[0], ADD_VALUES);
1520 } else {
1521 CHKERR VecSetValues(getFEMethod()->ts_F, row_data.getIndices().size(),
1522 &row_data.getIndices()[0], &nf[0], ADD_VALUES);
1523 }
1524
1525 } catch (const std::exception &ex) {
1526 std::ostringstream ss;
1527 ss << "throw in method: " << ex.what() << std::endl;
1528 SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
1529 }
1530
1532}
1533
1535 OpEshelbyDynamicMaterialMomentumLhs_dv(const std::string vel_field,
1536 const std::string field_name,
1537 BlockData &data,
1538 CommonData &common_data,
1539 Range *forcesonlyonentities_ptr)
1541 vel_field, field_name, data, common_data, forcesonlyonentities_ptr) {}
1542
1545 EntitiesFieldData::EntData &col_data, int gg) {
1547 int nb_col = col_data.getIndices().size();
1548 jac.clear();
1549 double *base_ptr = const_cast<double *>(&col_data.getN(gg)[0]);
1550 FTensor::Tensor0<double *> base(base_ptr, 1);
1551 double *diff_ptr =
1552 const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
1553 FTensor::Tensor1<double *, 3> diff(diff_ptr, &diff_ptr[1], &diff_ptr[2], 3);
1554 FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
1555 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1556 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1557 const int u = 3;
1559 &commonData.jacT[gg](0, u + 0), &commonData.jacT[gg](0, u + 1),
1560 &commonData.jacT[gg](0, u + 2), &commonData.jacT[gg](1, u + 0),
1561 &commonData.jacT[gg](1, u + 1), &commonData.jacT[gg](1, u + 2),
1562 &commonData.jacT[gg](2, u + 0), &commonData.jacT[gg](2, u + 1),
1563 &commonData.jacT[gg](2, u + 2));
1564 const int s = 3 + 3;
1566 // T* d000, T* d001, T* d002,
1567 // T* d010, T* d011, T* d012,
1568 // T* d020, T* d021, T* d022,
1569 // T* d100, T* d101, T* d102,
1570 // T* d110, T* d111, T* d112,
1571 // T* d120, T* d121, T* d122,
1572 // T* d200, T* d201, T* d202,
1573 // T* d210, T* d211, T* d212,
1574 // T* d220, T* d221, T* d222,
1575 &commonData.jacT[gg](0, s + 0), &commonData.jacT[gg](0, s + 1),
1576 &commonData.jacT[gg](0, s + 2), &commonData.jacT[gg](0, s + 3),
1577 &commonData.jacT[gg](0, s + 4), &commonData.jacT[gg](0, s + 5),
1578 &commonData.jacT[gg](0, s + 6), &commonData.jacT[gg](0, s + 7),
1579 &commonData.jacT[gg](0, s + 8), &commonData.jacT[gg](1, s + 0),
1580 &commonData.jacT[gg](1, s + 1), &commonData.jacT[gg](1, s + 2),
1581 &commonData.jacT[gg](1, s + 3), &commonData.jacT[gg](1, s + 4),
1582 &commonData.jacT[gg](1, s + 5), &commonData.jacT[gg](1, s + 6),
1583 &commonData.jacT[gg](1, s + 7), &commonData.jacT[gg](1, s + 8),
1584 &commonData.jacT[gg](2, s + 0), &commonData.jacT[gg](2, s + 1),
1585 &commonData.jacT[gg](2, s + 2), &commonData.jacT[gg](2, s + 3),
1586 &commonData.jacT[gg](2, s + 4), &commonData.jacT[gg](2, s + 5),
1587 &commonData.jacT[gg](2, s + 6), &commonData.jacT[gg](2, s + 7),
1588 &commonData.jacT[gg](2, s + 8));
1592 for (int dd = 0; dd < nb_col / 3; dd++) {
1593 t_jac(i, j) += t_mass1(i, j) * base * getFEMethod()->ts_a;
1594 t_jac(i, j) += t_mass3(i, j, k) * diff(k);
1595 ++base_ptr;
1596 ++diff_ptr;
1597 ++t_jac;
1598 }
1600}
1601
1603 OpEshelbyDynamicMaterialMomentumLhs_dx(const std::string vel_field,
1604 const std::string field_name,
1605 BlockData &data,
1606 CommonData &common_data,
1607 Range *forcesonlyonentities_ptr)
1609 vel_field, field_name, data, common_data, forcesonlyonentities_ptr) {}
1610
1613 EntitiesFieldData::EntData &col_data, int gg) {
1615 int nb_col = col_data.getIndices().size();
1616 jac.clear();
1617 double *diff_ptr =
1618 const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
1619 FTensor::Tensor1<double *, 3> diff(diff_ptr, &diff_ptr[1], &diff_ptr[2], 3);
1620 FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
1621 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1622 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1623 const int s = 3 + 3 + 9;
1625 // T* d000, T* d001, T* d002,
1626 // T* d010, T* d011, T* d012,
1627 // T* d020, T* d021, T* d022,
1628 // T* d100, T* d101, T* d102,
1629 // T* d110, T* d111, T* d112,
1630 // T* d120, T* d121, T* d122,
1631 // T* d200, T* d201, T* d202,
1632 // T* d210, T* d211, T* d212,
1633 // T* d220, T* d221, T* d222,
1634 &commonData.jacT[gg](0, s + 0), &commonData.jacT[gg](0, s + 1),
1635 &commonData.jacT[gg](0, s + 2), &commonData.jacT[gg](0, s + 3),
1636 &commonData.jacT[gg](0, s + 4), &commonData.jacT[gg](0, s + 5),
1637 &commonData.jacT[gg](0, s + 6), &commonData.jacT[gg](0, s + 7),
1638 &commonData.jacT[gg](0, s + 8), &commonData.jacT[gg](1, s + 0),
1639 &commonData.jacT[gg](1, s + 1), &commonData.jacT[gg](1, s + 2),
1640 &commonData.jacT[gg](1, s + 3), &commonData.jacT[gg](1, s + 4),
1641 &commonData.jacT[gg](1, s + 5), &commonData.jacT[gg](1, s + 6),
1642 &commonData.jacT[gg](1, s + 7), &commonData.jacT[gg](1, s + 8),
1643 &commonData.jacT[gg](2, s + 0), &commonData.jacT[gg](2, s + 1),
1644 &commonData.jacT[gg](2, s + 2), &commonData.jacT[gg](2, s + 3),
1645 &commonData.jacT[gg](2, s + 4), &commonData.jacT[gg](2, s + 5),
1646 &commonData.jacT[gg](2, s + 6), &commonData.jacT[gg](2, s + 7),
1647 &commonData.jacT[gg](2, s + 8));
1651 for (int dd = 0; dd < nb_col / 3; dd++) {
1652 t_jac(i, j) += t_mass3(i, j, k) * diff(k);
1653 ++diff_ptr;
1654 ++t_jac;
1655 }
1657}
1658
1660 OpEshelbyDynamicMaterialMomentumLhs_dX(const std::string vel_field,
1661 const std::string field_name,
1662 BlockData &data,
1663 CommonData &common_data,
1664 Range *forcesonlyonentities_ptr)
1666 vel_field, field_name, data, common_data, forcesonlyonentities_ptr) {}
1667
1670 EntitiesFieldData::EntData &col_data, int gg) {
1672 int nb_col = col_data.getIndices().size();
1673 jac.clear();
1674 double *diff_ptr =
1675 const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
1676 FTensor::Tensor1<double *, 3> diff(diff_ptr, &diff_ptr[1], &diff_ptr[2], 3);
1677 FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
1678 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1679 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1680 const int s = 3 + 3 + 9 + 9;
1682 // T* d000, T* d001, T* d002,
1683 // T* d010, T* d011, T* d012,
1684 // T* d020, T* d021, T* d022,
1685 // T* d100, T* d101, T* d102,
1686 // T* d110, T* d111, T* d112,
1687 // T* d120, T* d121, T* d122,
1688 // T* d200, T* d201, T* d202,
1689 // T* d210, T* d211, T* d212,
1690 // T* d220, T* d221, T* d222,
1691 &commonData.jacT[gg](0, s + 0), &commonData.jacT[gg](0, s + 1),
1692 &commonData.jacT[gg](0, s + 2), &commonData.jacT[gg](0, s + 3),
1693 &commonData.jacT[gg](0, s + 4), &commonData.jacT[gg](0, s + 5),
1694 &commonData.jacT[gg](0, s + 6), &commonData.jacT[gg](0, s + 7),
1695 &commonData.jacT[gg](0, s + 8), &commonData.jacT[gg](1, s + 0),
1696 &commonData.jacT[gg](1, s + 1), &commonData.jacT[gg](1, s + 2),
1697 &commonData.jacT[gg](1, s + 3), &commonData.jacT[gg](1, s + 4),
1698 &commonData.jacT[gg](1, s + 5), &commonData.jacT[gg](1, s + 6),
1699 &commonData.jacT[gg](1, s + 7), &commonData.jacT[gg](1, s + 8),
1700 &commonData.jacT[gg](2, s + 0), &commonData.jacT[gg](2, s + 1),
1701 &commonData.jacT[gg](2, s + 2), &commonData.jacT[gg](2, s + 3),
1702 &commonData.jacT[gg](2, s + 4), &commonData.jacT[gg](2, s + 5),
1703 &commonData.jacT[gg](2, s + 6), &commonData.jacT[gg](2, s + 7),
1704 &commonData.jacT[gg](2, s + 8));
1708 for (int dd = 0; dd < nb_col / 3; dd++) {
1709 t_jac(i, j) += t_mass3(i, j, k) * diff(k);
1710 ++diff_ptr;
1711 ++t_jac;
1712 }
1714}
1715
1717 MoFEM::Interface &m_field, TS _ts, const std::string velocity_field,
1718 const std::string spatial_position_field)
1719 : mField(m_field), tS(_ts), velocityField(velocity_field),
1720 spatialPositionField(spatial_position_field), jacobianLag(-1) {}
1721
1724
1725 switch (ts_ctx) {
1726 case CTX_TSSETIFUNCTION: {
1727 snes_f = ts_F;
1728 // FIXME: This global scattering because Kuu problem and Dynamic problem
1729 // not share partitions. Both problem should use the same partitioning to
1730 // resolve this problem.
1731 CHKERR mField.getInterface<VecManager>()->setGlobalGhostVector(
1732 problemPtr, COL, ts_u, INSERT_VALUES, SCATTER_REVERSE);
1733 CHKERR mField.getInterface<VecManager>()->setOtherGlobalGhostVector(
1734 problemPtr, velocityField, "DOT_" + velocityField, COL, ts_u_t,
1735 INSERT_VALUES, SCATTER_REVERSE);
1736 CHKERR mField.getInterface<VecManager>()->setOtherGlobalGhostVector(
1737 problemPtr, spatialPositionField, "DOT_" + spatialPositionField, COL,
1738 ts_u_t, INSERT_VALUES, SCATTER_REVERSE);
1739 break;
1740 }
1741 case CTX_TSSETIJACOBIAN: {
1742 snes_B = ts_B;
1743 break;
1744 }
1745 default:
1746 break;
1747 }
1748
1750}
1751
1754 //
1755 // SNES snes;
1756 // CHKERR TSGetSNES(tS,&snes);
1757 // CHKERR SNESSetLagJacobian(snes,jacobianLag);
1759}
1760
1763
1764 Range added_tets;
1766 mField, BLOCKSET | BODYFORCESSET, it)) {
1767 int id = it->getMeshsetId();
1768 EntityHandle meshset = it->getMeshset();
1769 CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET,
1770 setOfBlocks[id].tEts, true);
1771 added_tets.merge(setOfBlocks[id].tEts);
1772 Block_BodyForces mydata;
1773 CHKERR it->getAttributeDataStructure(mydata);
1774 setOfBlocks[id].rho0 = mydata.data.density;
1775 setOfBlocks[id].a0.resize(3);
1776 setOfBlocks[id].a0[0] = mydata.data.acceleration_x;
1777 setOfBlocks[id].a0[1] = mydata.data.acceleration_y;
1778 setOfBlocks[id].a0[2] = mydata.data.acceleration_z;
1779 // std::cerr << setOfBlocks[id].tEts << std::endl;
1780 }
1781
1783 mField, BLOCKSET | MAT_ELASTICSET, it)) {
1784 Mat_Elastic mydata;
1785 CHKERR it->getAttributeDataStructure(mydata);
1786 if (mydata.data.User1 == 0)
1787 continue;
1788 Range tets;
1789 EntityHandle meshset = it->getMeshset();
1790 CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets, true);
1791 tets = subtract(tets, added_tets);
1792 if (tets.empty())
1793 continue;
1794 int id = it->getMeshsetId();
1795 setOfBlocks[-id].tEts = tets;
1796 setOfBlocks[-id].rho0 = mydata.data.User1;
1797 setOfBlocks[-id].a0.resize(3);
1798 setOfBlocks[-id].a0[0] = mydata.data.User2;
1799 setOfBlocks[-id].a0[1] = mydata.data.User3;
1800 setOfBlocks[-id].a0[2] = mydata.data.User4;
1801 // std::cerr << setOfBlocks[id].tEts << std::endl;
1802 }
1803
1805}
1806
1808 MoFEM::Interface &m_field,
1809 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr) {
1811
1812 if (!block_sets_ptr)
1813 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1814 "Pointer to block of sets is null");
1815
1817 m_field, BLOCKSET | BODYFORCESSET, it)) {
1818 Block_BodyForces mydata;
1819 CHKERR it->getAttributeDataStructure(mydata);
1820 int id = it->getMeshsetId();
1821 auto &block_data = (*block_sets_ptr)[id];
1822 EntityHandle meshset = it->getMeshset();
1823 CHKERR m_field.get_moab().get_entities_by_dimension(meshset, 3,
1824 block_data.tEts, true);
1825 block_data.rho0 = mydata.data.density;
1826 block_data.a0.resize(3);
1827 block_data.a0[0] = mydata.data.acceleration_x;
1828 block_data.a0[1] = mydata.data.acceleration_y;
1829 block_data.a0[2] = mydata.data.acceleration_z;
1830 }
1831
1833}
1834
1836 string element_name, string velocity_field_name,
1837 string spatial_position_field_name, string material_position_field_name,
1838 bool ale, BitRefLevel bit) {
1840
1841 //
1842
1843 CHKERR mField.add_finite_element(element_name, MF_ZERO);
1845 velocity_field_name);
1847 velocity_field_name);
1849 velocity_field_name);
1851 element_name, spatial_position_field_name);
1853 element_name, spatial_position_field_name);
1855 element_name, spatial_position_field_name);
1856 if (mField.check_field(material_position_field_name)) {
1857 if (ale) {
1859 element_name, material_position_field_name);
1861 element_name, material_position_field_name);
1863 element_name, "DOT_" + material_position_field_name);
1864 }
1866 element_name, material_position_field_name);
1867 }
1869 element_name, "DOT_" + velocity_field_name);
1871 element_name, "DOT_" + spatial_position_field_name);
1872
1873 Range tets;
1874 if (bit.any()) {
1875 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1876 bit, BitRefLevel().set(), MBTET, tets);
1877 }
1878
1879 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
1880 for (; sit != setOfBlocks.end(); sit++) {
1881 Range add_tets = sit->second.tEts;
1882 if (!tets.empty()) {
1883 add_tets = intersect(add_tets, tets);
1884 }
1886 element_name);
1887 }
1888
1890}
1891
1893 string element_name, string velocity_field_name,
1894 string spatial_position_field_name, string material_position_field_name,
1895 bool ale, BitRefLevel bit) {
1897
1898 //
1899
1900 CHKERR mField.add_finite_element(element_name, MF_ZERO);
1902 velocity_field_name);
1904 velocity_field_name);
1906 velocity_field_name);
1908 element_name, spatial_position_field_name);
1910 element_name, spatial_position_field_name);
1911 if (mField.check_field(material_position_field_name)) {
1912 if (ale) {
1914 element_name, material_position_field_name);
1916 element_name, "DOT_" + material_position_field_name);
1917 }
1919 element_name, material_position_field_name);
1920 }
1922 element_name, "DOT_" + velocity_field_name);
1924 element_name, "DOT_" + spatial_position_field_name);
1925
1926 Range tets;
1927 if (bit.any()) {
1928 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1929 bit, BitRefLevel().set(), MBTET, tets);
1930 }
1931
1932 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
1933 for (; sit != setOfBlocks.end(); sit++) {
1934 Range add_tets = sit->second.tEts;
1935 if (!tets.empty()) {
1936 add_tets = intersect(add_tets, tets);
1937 }
1939 element_name);
1940 }
1941
1943}
1944
1946 string element_name, string velocity_field_name,
1947 string spatial_position_field_name, string material_position_field_name,
1948 bool ale, BitRefLevel bit, Range *intersected) {
1950
1951 //
1952
1953 CHKERR mField.add_finite_element(element_name, MF_ZERO);
1955 velocity_field_name);
1957 velocity_field_name);
1959 element_name, spatial_position_field_name);
1961 element_name, spatial_position_field_name);
1962 if (mField.check_field(material_position_field_name)) {
1963 if (ale) {
1965 element_name, material_position_field_name);
1967 element_name, material_position_field_name);
1969 element_name, "DOT_" + material_position_field_name);
1970 }
1972 element_name, material_position_field_name);
1973 }
1975 element_name, "DOT_" + velocity_field_name);
1977 element_name, "DOT_" + spatial_position_field_name);
1978
1979 Range tets;
1980 if (bit.any()) {
1981 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1982 bit, BitRefLevel().set(), MBTET, tets);
1983 }
1984 if (intersected != NULL) {
1985 if (tets.empty()) {
1986 tets = *intersected;
1987 } else {
1988 tets = intersect(*intersected, tets);
1989 }
1990 }
1991
1992 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
1993 for (; sit != setOfBlocks.end(); sit++) {
1994 Range add_tets = sit->second.tEts;
1995 if (!tets.empty()) {
1996 add_tets = intersect(add_tets, tets);
1997 }
1999 element_name);
2000 }
2001
2003}
2004
2006 string velocity_field_name, string spatial_position_field_name,
2007 string material_position_field_name, bool ale, bool linear) {
2009
2010 commonData.spatialPositions = spatial_position_field_name;
2011 commonData.meshPositions = material_position_field_name;
2012 commonData.spatialVelocities = velocity_field_name;
2013 commonData.lInear = linear;
2014
2015 // Rhs
2016 feMassRhs.getOpPtrVector().push_back(
2017 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2018 feMassRhs.getOpPtrVector().push_back(
2019 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2020 feMassRhs.getOpPtrVector().push_back(
2021 new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2022 feMassRhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2023 "DOT_" + spatial_position_field_name, commonData));
2024 if (mField.check_field(material_position_field_name)) {
2025 feMassRhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2026 material_position_field_name, commonData));
2027 if (ale) {
2028 feMassRhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2029 "DOT_" + material_position_field_name, commonData));
2030 } else {
2031 feMassRhs.meshPositionsFieldName = material_position_field_name;
2032 }
2033 }
2034 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
2035 for (; sit != setOfBlocks.end(); sit++) {
2036 feMassRhs.getOpPtrVector().push_back(
2037 new OpMassJacobian(spatial_position_field_name, sit->second, commonData,
2038 methodsOp, tAg, false));
2039 feMassRhs.getOpPtrVector().push_back(
2040 new OpMassRhs(spatial_position_field_name, sit->second, commonData));
2041 }
2042
2043 // Lhs
2044 feMassLhs.getOpPtrVector().push_back(
2045 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2046 feMassLhs.getOpPtrVector().push_back(
2047 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2048 feMassLhs.getOpPtrVector().push_back(
2049 new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2050 feMassLhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2051 "DOT_" + spatial_position_field_name, commonData));
2052 if (mField.check_field(material_position_field_name)) {
2053 feMassLhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2054 material_position_field_name, commonData));
2055 if (ale) {
2056 feMassLhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2057 "DOT_" + material_position_field_name, commonData));
2058 } else {
2059 feMassLhs.meshPositionsFieldName = material_position_field_name;
2060 }
2061 }
2062 sit = setOfBlocks.begin();
2063 for (; sit != setOfBlocks.end(); sit++) {
2064 feMassLhs.getOpPtrVector().push_back(
2065 new OpMassJacobian(spatial_position_field_name, sit->second, commonData,
2066 methodsOp, tAg, true));
2067 feMassLhs.getOpPtrVector().push_back(
2068 new OpMassLhs_dM_dv(spatial_position_field_name, velocity_field_name,
2069 sit->second, commonData));
2070 feMassLhs.getOpPtrVector().push_back(new OpMassLhs_dM_dx(
2071 spatial_position_field_name, spatial_position_field_name, sit->second,
2072 commonData));
2073 if (mField.check_field(material_position_field_name)) {
2074 if (ale) {
2075 feMassLhs.getOpPtrVector().push_back(new OpMassLhs_dM_dX(
2076 spatial_position_field_name, material_position_field_name,
2077 sit->second, commonData));
2078 } else {
2079 feMassLhs.meshPositionsFieldName = material_position_field_name;
2080 }
2081 }
2082 }
2083
2084 // Energy
2085 feEnergy.getOpPtrVector().push_back(
2086 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2087 feEnergy.getOpPtrVector().push_back(
2088 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2089 if (mField.check_field(material_position_field_name)) {
2090 feEnergy.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2091 material_position_field_name, commonData));
2092 feEnergy.meshPositionsFieldName = material_position_field_name;
2093 }
2094 sit = setOfBlocks.begin();
2095 for (; sit != setOfBlocks.end(); sit++) {
2096 feEnergy.getOpPtrVector().push_back(new OpEnergy(
2097 spatial_position_field_name, sit->second, commonData, feEnergy.V));
2098 }
2099
2101}
2102
2104 string velocity_field_name, string spatial_position_field_name,
2105 string material_position_field_name, bool ale) {
2107
2108 commonData.spatialPositions = spatial_position_field_name;
2109 commonData.meshPositions = material_position_field_name;
2110 commonData.spatialVelocities = velocity_field_name;
2111
2112 // Rhs
2113 feVelRhs.getOpPtrVector().push_back(
2114 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2115 feVelRhs.getOpPtrVector().push_back(
2116 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2117 feVelRhs.getOpPtrVector().push_back(
2118 new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2119 if (mField.check_field(material_position_field_name)) {
2120 feVelRhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2121 "DOT_" + spatial_position_field_name, commonData));
2122 feVelRhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2123 material_position_field_name, commonData));
2124 if (ale) {
2125 feVelRhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2126 "DOT_" + material_position_field_name, commonData));
2127 } else {
2128 feVelRhs.meshPositionsFieldName = material_position_field_name;
2129 }
2130 }
2131 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
2132 for (; sit != setOfBlocks.end(); sit++) {
2133 feVelRhs.getOpPtrVector().push_back(new OpVelocityJacobian(
2134 velocity_field_name, sit->second, commonData, tAg, false));
2135 feVelRhs.getOpPtrVector().push_back(
2136 new OpVelocityRhs(velocity_field_name, sit->second, commonData));
2137 }
2138
2139 // Lhs
2140 feVelLhs.getOpPtrVector().push_back(
2141 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2142 feVelLhs.getOpPtrVector().push_back(
2143 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2144 feVelLhs.getOpPtrVector().push_back(
2145 new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2146 if (mField.check_field(material_position_field_name)) {
2147 feVelLhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2148 "DOT_" + spatial_position_field_name, commonData));
2149 feVelLhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2150 material_position_field_name, commonData));
2151 if (ale) {
2152 feVelLhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2153 "DOT_" + material_position_field_name, commonData));
2154 } else {
2155 feVelLhs.meshPositionsFieldName = material_position_field_name;
2156 }
2157 }
2158 sit = setOfBlocks.begin();
2159 for (; sit != setOfBlocks.end(); sit++) {
2160 feVelLhs.getOpPtrVector().push_back(new OpVelocityJacobian(
2161 velocity_field_name, sit->second, commonData, tAg));
2162 feVelLhs.getOpPtrVector().push_back(new OpVelocityLhs_dV_dv(
2163 velocity_field_name, velocity_field_name, sit->second, commonData));
2164 feVelLhs.getOpPtrVector().push_back(new OpVelocityLhs_dV_dx(
2165 velocity_field_name, spatial_position_field_name, sit->second,
2166 commonData));
2167 if (mField.check_field(material_position_field_name)) {
2168 if (ale) {
2169 feVelLhs.getOpPtrVector().push_back(new OpVelocityLhs_dV_dX(
2170 velocity_field_name, material_position_field_name, sit->second,
2171 commonData));
2172 } else {
2173 feVelLhs.meshPositionsFieldName = material_position_field_name;
2174 }
2175 }
2176 }
2177
2179}
2180
2182 string velocity_field_name, string spatial_position_field_name,
2183 string material_position_field_name, Range *forces_on_entities_ptr) {
2185
2186 commonData.spatialPositions = spatial_position_field_name;
2187 commonData.meshPositions = material_position_field_name;
2188 commonData.spatialVelocities = velocity_field_name;
2189
2190 // Rhs
2191 feTRhs.getOpPtrVector().push_back(
2192 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2193 feTRhs.getOpPtrVector().push_back(
2194 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2195 feTRhs.getOpPtrVector().push_back(
2196 new OpGetCommonDataAtGaussPts(material_position_field_name, commonData));
2197 feTRhs.getOpPtrVector().push_back(
2198 new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2199
2200 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
2201 for (; sit != setOfBlocks.end(); sit++) {
2202 feTRhs.getOpPtrVector().push_back(
2204 material_position_field_name, sit->second, commonData, tAg, false));
2205 feTRhs.getOpPtrVector().push_back(new OpEshelbyDynamicMaterialMomentumRhs(
2206 material_position_field_name, sit->second, commonData,
2207 forces_on_entities_ptr));
2208 }
2209
2210 // Lhs
2211 feTLhs.getOpPtrVector().push_back(
2212 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2213 feTLhs.getOpPtrVector().push_back(
2214 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2215 feTLhs.getOpPtrVector().push_back(
2216 new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2217 if (mField.check_field(material_position_field_name)) {
2218 feTLhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2219 material_position_field_name, commonData));
2220 }
2221 sit = setOfBlocks.begin();
2222 for (; sit != setOfBlocks.end(); sit++) {
2223 feTLhs.getOpPtrVector().push_back(
2225 material_position_field_name, sit->second, commonData, tAg));
2226 feTLhs.getOpPtrVector().push_back(
2228 material_position_field_name, velocity_field_name, sit->second,
2229 commonData, forces_on_entities_ptr));
2230 feTLhs.getOpPtrVector().push_back(
2232 material_position_field_name, spatial_position_field_name,
2233 sit->second, commonData, forces_on_entities_ptr));
2234 feTLhs.getOpPtrVector().push_back(
2236 material_position_field_name, material_position_field_name,
2237 sit->second, commonData, forces_on_entities_ptr));
2238 }
2239
2241}
2242
2244 string velocity_field_name, string spatial_position_field_name,
2245 string material_position_field_name, bool linear) {
2247
2248 commonData.spatialPositions = spatial_position_field_name;
2249 commonData.meshPositions = material_position_field_name;
2250 commonData.spatialVelocities = velocity_field_name;
2251 commonData.lInear = linear;
2252
2253 // Rhs
2254 feMassRhs.getOpPtrVector().push_back(
2255 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2256 feMassRhs.getOpPtrVector().push_back(
2257 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2258 feMassRhs.getOpPtrVector().push_back(
2259 new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2260 if (mField.check_field(material_position_field_name)) {
2261 feMassRhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2262 material_position_field_name, commonData));
2263 feMassRhs.meshPositionsFieldName = material_position_field_name;
2264 }
2265 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
2266 for (; sit != setOfBlocks.end(); sit++) {
2267 feMassRhs.getOpPtrVector().push_back(
2268 new OpMassJacobian(spatial_position_field_name, sit->second, commonData,
2269 methodsOp, tAg, false));
2270 feMassRhs.getOpPtrVector().push_back(
2271 new OpMassRhs(spatial_position_field_name, sit->second, commonData));
2272 }
2273
2274 // Lhs
2275 feMassLhs.getOpPtrVector().push_back(
2276 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2277 feMassLhs.getOpPtrVector().push_back(
2278 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2279 feMassLhs.getOpPtrVector().push_back(
2280 new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2281 if (mField.check_field(material_position_field_name)) {
2282 feMassLhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2283 material_position_field_name, commonData));
2284 feMassLhs.meshPositionsFieldName = material_position_field_name;
2285 }
2286 sit = setOfBlocks.begin();
2287 for (; sit != setOfBlocks.end(); sit++) {
2288 feMassLhs.getOpPtrVector().push_back(
2289 new OpMassJacobian(spatial_position_field_name, sit->second, commonData,
2290 methodsOp, tAg, true));
2291 feMassLhs.getOpPtrVector().push_back(new OpMassLhs_dM_dv(
2292 spatial_position_field_name, spatial_position_field_name, sit->second,
2293 commonData));
2294 if (mField.check_field(material_position_field_name)) {
2295 feMassLhs.meshPositionsFieldName = material_position_field_name;
2296 }
2297 }
2298
2299 // Aux Lhs
2300 feMassAuxLhs.getOpPtrVector().push_back(
2301 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2302 feMassAuxLhs.getOpPtrVector().push_back(
2303 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2304 feMassAuxLhs.getOpPtrVector().push_back(
2305 new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2306 if (mField.check_field(material_position_field_name)) {
2307 feMassAuxLhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2308 material_position_field_name, commonData));
2309 feMassAuxLhs.meshPositionsFieldName = material_position_field_name;
2310 }
2311 sit = setOfBlocks.begin();
2312 for (; sit != setOfBlocks.end(); sit++) {
2313 feMassAuxLhs.getOpPtrVector().push_back(
2314 new OpMassJacobian(spatial_position_field_name, sit->second, commonData,
2315 methodsOp, tAg, true));
2316 feMassAuxLhs.getOpPtrVector().push_back(new OpMassLhs_dM_dx(
2317 spatial_position_field_name, spatial_position_field_name, sit->second,
2318 commonData));
2319 if (mField.check_field(material_position_field_name)) {
2320 feMassAuxLhs.meshPositionsFieldName = material_position_field_name;
2321 }
2322 }
2323
2324 // Energy E=0.5*rho*v*v
2325 feEnergy.getOpPtrVector().push_back(
2326 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2327 feEnergy.getOpPtrVector().push_back(
2328 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2329 if (mField.check_field(material_position_field_name)) {
2330 feEnergy.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2331 material_position_field_name, commonData));
2332 feEnergy.meshPositionsFieldName = material_position_field_name;
2333 }
2334 sit = setOfBlocks.begin();
2335 for (; sit != setOfBlocks.end(); sit++) {
2336 feEnergy.getOpPtrVector().push_back(new OpEnergy(
2337 spatial_position_field_name, sit->second, commonData, feEnergy.V));
2338 }
2339
2341}
2342
2345 if (iNitialized) {
2346
2347 CHKERR dEstroy();
2348 CHKERRABORT(PETSC_COMM_WORLD, ierr);
2349 }
2350}
2351
2354 if (!iNitialized) {
2355
2356#if PETSC_VERSION_GE(3, 5, 3)
2357 CHKERR MatCreateVecs(K, &u, &Ku);
2358 CHKERR MatCreateVecs(M, &v, &Mv);
2359#else
2360 CHKERR MatGetVecs(K, &u, &Ku);
2361 CHKERR MatGetVecs(M, &v, &Mv);
2362#endif
2363 CHKERR MatDuplicate(K, MAT_SHARE_NONZERO_PATTERN, &barK);
2364 iNitialized = true;
2365 }
2367}
2368
2371 if (iNitialized) {
2372
2373 CHKERR VecDestroy(&u);
2374 CHKERR VecDestroy(&Ku);
2375 CHKERR VecDestroy(&v);
2376 CHKERR VecDestroy(&Mv);
2377 CHKERR MatDestroy(&barK);
2378 iNitialized = false;
2379 }
2381}
2382
2385
2386 if (!initPC) {
2387 MPI_Comm comm;
2388 CHKERR PetscObjectGetComm((PetscObject)shellMat, &comm);
2389 CHKERR PCCreate(comm, &pC);
2390 initPC = true;
2391 }
2393}
2394
2397
2398 if (initPC) {
2399 CHKERR PCDestroy(&pC);
2400 initPC = false;
2401 }
2403}
2404
2406 MoFEM::Interface &m_field)
2407 : mField(m_field) {}
2408
2411
2412 if (ts_ctx != CTX_TSSETIFUNCTION) {
2413 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
2414 "It is used to residual of velocities");
2415 }
2416 if (!shellMatCtx->iNitialized) {
2417 CHKERR shellMatCtx->iNit();
2418 }
2419 // Note velocities calculate from displacements are stroed in shellMatCtx->u
2420 CHKERR VecScatterBegin(shellMatCtx->scatterU, ts_u_t, shellMatCtx->u,
2421 INSERT_VALUES, SCATTER_FORWARD);
2422 CHKERR VecScatterEnd(shellMatCtx->scatterU, ts_u_t, shellMatCtx->u,
2423 INSERT_VALUES, SCATTER_FORWARD);
2424 CHKERR VecScatterBegin(shellMatCtx->scatterV, ts_u, shellMatCtx->v,
2425 INSERT_VALUES, SCATTER_FORWARD);
2426 CHKERR VecScatterEnd(shellMatCtx->scatterV, ts_u, shellMatCtx->v,
2427 INSERT_VALUES, SCATTER_FORWARD);
2428 CHKERR VecAXPY(shellMatCtx->v, -1, shellMatCtx->u);
2429 CHKERR VecScatterBegin(shellMatCtx->scatterV, shellMatCtx->v, ts_F,
2430 ADD_VALUES, SCATTER_REVERSE);
2431 CHKERR VecScatterEnd(shellMatCtx->scatterV, shellMatCtx->v, ts_F, ADD_VALUES,
2432 SCATTER_REVERSE);
2433 // VecView(shellMatCtx->v,PETSC_VIEWER_STDOUT_WORLD);
2434
2436}
2437
2438#ifdef __DIRICHLET_HPP__
2439
2440ConvectiveMassElement::ShellMatrixElement::ShellMatrixElement(
2441 MoFEM::Interface &m_field)
2442 : mField(m_field) {}
2443
2444MoFEMErrorCode ConvectiveMassElement::ShellMatrixElement::preProcess() {
2446
2447 if (ts_ctx != CTX_TSSETIJACOBIAN) {
2448 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
2449 "It is used to calculate shell matrix only");
2450 }
2451
2452 shellMatCtx->ts_a = ts_a;
2453 DirichletBcPtr->copyTs(*((TSMethod *)this)); // copy context for TSMethod
2454
2455 DirichletBcPtr->dIag = 1;
2456 DirichletBcPtr->ts_B = shellMatCtx->K;
2457 CHKERR MatZeroEntries(shellMatCtx->K);
2458 CHKERR mField.problem_basic_method_preProcess(problemName, *DirichletBcPtr);
2459 LoopsToDoType::iterator itk = loopK.begin();
2460 for (; itk != loopK.end(); itk++) {
2461 itk->second->copyTs(*((TSMethod *)this));
2462 itk->second->ts_B = shellMatCtx->K;
2463 CHKERR mField.loop_finite_elements(problemName, itk->first, *itk->second);
2464 }
2465 LoopsToDoType::iterator itam = loopAuxM.begin();
2466 for (; itam != loopAuxM.end(); itam++) {
2467 itam->second->copyTs(*((TSMethod *)this));
2468 itam->second->ts_B = shellMatCtx->K;
2469 CHKERR mField.loop_finite_elements(problemName, itam->first, *itam->second);
2470 }
2471 CHKERR mField.problem_basic_method_postProcess(problemName, *DirichletBcPtr);
2472 CHKERR MatAssemblyBegin(shellMatCtx->K, MAT_FINAL_ASSEMBLY);
2473 CHKERR MatAssemblyEnd(shellMatCtx->K, MAT_FINAL_ASSEMBLY);
2474
2475 DirichletBcPtr->dIag = 0;
2476 DirichletBcPtr->ts_B = shellMatCtx->M;
2477 CHKERR MatZeroEntries(shellMatCtx->M);
2478 // CHKERR mField.problem_basic_method_preProcess(problemName,*DirichletBcPtr);
2479 LoopsToDoType::iterator itm = loopM.begin();
2480 for (; itm != loopM.end(); itm++) {
2481 itm->second->copyTs(*((TSMethod *)this));
2482 itm->second->ts_B = shellMatCtx->M;
2483 CHKERR mField.loop_finite_elements(problemName, itm->first, *itm->second);
2484 }
2485 CHKERR mField.problem_basic_method_postProcess(problemName, *DirichletBcPtr);
2486 CHKERR MatAssemblyBegin(shellMatCtx->M, MAT_FINAL_ASSEMBLY);
2487 CHKERR MatAssemblyEnd(shellMatCtx->M, MAT_FINAL_ASSEMBLY);
2488
2489 // barK
2490 CHKERR MatZeroEntries(shellMatCtx->barK);
2491 CHKERR MatCopy(shellMatCtx->K, shellMatCtx->barK, SAME_NONZERO_PATTERN);
2492 CHKERR MatAXPY(shellMatCtx->barK, ts_a, shellMatCtx->M, SAME_NONZERO_PATTERN);
2493 CHKERR MatAssemblyBegin(shellMatCtx->barK, MAT_FINAL_ASSEMBLY);
2494 CHKERR MatAssemblyEnd(shellMatCtx->barK, MAT_FINAL_ASSEMBLY);
2495
2496 // Matrix View
2497 // MatView(shellMatCtx->barK,PETSC_VIEWER_DRAW_WORLD);//PETSC_VIEWER_STDOUT_WORLD);
2498 // std::string wait;
2499 // std::cin >> wait;
2500
2502}
2503
2504#endif //__DIRICHLET_HPP__
static Index< 'M', 3 > M
Operators and data structures for mass and convective mass element.
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr double a
@ COL
Definition: definitions.h:123
@ MF_ZERO
Definition: definitions.h:98
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ BODYFORCESSET
block name is "BODY_FORCES"
Definition: definitions.h:162
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
Definition: definitions.h:159
@ BLOCKSET
Definition: definitions.h:148
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
constexpr int order
@ F
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual bool check_field(const std::string &name) const =0
check if field is in database
virtual MoFEMErrorCode problem_basic_method_postProcess(const Problem *problem_ptr, BasicMethod &method, int verb=DEFAULT_VERBOSITY)=0
Set data for BasicMethod.
virtual MoFEMErrorCode loop_finite_elements(const std::string problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
auto bit
set bit
constexpr double a0
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
MoFEM::TsCtx * ts_ctx
Definition: level_set.cpp:1884
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
UBlasVector< int > VectorInt
Definition: Types.hpp:67
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
Definition: Templates.hpp:1608
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
auto getFTensor2FromArray3by3(ublas::matrix< T, L, A > &data, const FTensor::Number< S > &, const size_t rr, const size_t cc=0)
Definition: Templates.hpp:1359
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1557
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
constexpr AssemblyType A
double h
double rho
Definition: plastic.cpp:191
double H
Hardening.
Definition: plastic.cpp:175
constexpr auto field_name
constexpr double g
data for calculation inertia forces
common data used by volume elements
std::vector< std::vector< double * > > jacTRowPtr
std::map< std::string, std::vector< VectorDouble > > dataAtGaussPts
std::vector< std::vector< double * > > jacVelRowPtr
std::map< std::string, std::vector< MatrixDouble > > gradAtGaussPts
std::vector< std::vector< double * > > jacMassRowPtr
int getRule(int order)
it is used to calculate nb. of Gauss integration points
MyVolumeFE(MoFEM::Interface &m_field)
OpEnergy(const std::string field_name, BlockData &data, CommonData &common_data, SmartPetscObj< Vec > v)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpEshelbyDynamicMaterialMomentumJacobian(const std::string field_name, BlockData &data, CommonData &common_data, int tag, bool jacobian=true)
OpEshelbyDynamicMaterialMomentumLhs_dX(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data, Range *forcesonlyonentities_ptr)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
OpEshelbyDynamicMaterialMomentumLhs_dv(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data, Range *forcesonlyonentities_ptr)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
OpEshelbyDynamicMaterialMomentumLhs_dx(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data, Range *forcesonlyonentities_ptr)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpEshelbyDynamicMaterialMomentumRhs(const std::string field_name, BlockData &data, CommonData &common_data, Range *forcesonlyonentities_ptr)
OpGetCommonDataAtGaussPts(const std::string field_name, CommonData &common_data)
OpGetDataAtGaussPts(const std::string field_name, std::vector< VectorDouble > &values_at_gauss_pts, std::vector< MatrixDouble > &gardient_at_gauss_pts)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
operator calculating deformation gradient
OpMassJacobian(const std::string field_name, BlockData &data, CommonData &common_data, boost::ptr_vector< MethodForForceScaling > &methods_op, int tag, bool linear=false)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpMassLhs_dM_dX(const std::string field_name, const std::string col_field, BlockData &data, CommonData &common_data)
MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
OpMassLhs_dM_dv(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data, Range *forcesonlyonentities_ptr=NULL)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
OpMassLhs_dM_dx(const std::string field_name, const std::string col_field, BlockData &data, CommonData &common_data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpMassRhs(const std::string field_name, BlockData &data, CommonData &common_data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpVelocityJacobian(const std::string field_name, BlockData &data, CommonData &common_data, int tag, bool jacobian=true)
OpVelocityLhs_dV_dX(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
OpVelocityLhs_dV_dv(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
OpVelocityLhs_dV_dx(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
OpVelocityRhs(const std::string field_name, BlockData &data, CommonData &common_data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
MoFEMErrorCode preProcess()
Calculate inconsistency between approximation of velocities and velocities calculated from displaceme...
UpdateAndControl(MoFEM::Interface &m_field, TS _ts, const std::string velocity_field, const std::string spatial_position_field)
MoFEMErrorCode preProcess()
Scatter values from t_u_dt on the fields.
structure grouping operators and data used for calculation of mass (convective) element \ nonlinear_e...
ConvectiveMassElement(MoFEM::Interface &m_field, short int tag)
MoFEMErrorCode setVelocityOperators(string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool ale=false)
MyVolumeFE feEnergy
calculate kinetic energy
MoFEMErrorCode setShellMatrixMassOperators(string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool linear=false)
MyVolumeFE feVelRhs
calculate right hand side for tetrahedral elements
MoFEMErrorCode addVelocityElement(string element_name, string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool ale=false, BitRefLevel bit=BitRefLevel())
MoFEMErrorCode addConvectiveMassElement(string element_name, string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool ale=false, BitRefLevel bit=BitRefLevel())
MoFEMErrorCode setKinematicEshelbyOperators(string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", Range *forces_on_entities_ptr=NULL)
MoFEMErrorCode addEshelbyDynamicMaterialMomentum(string element_name, string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool ale=false, BitRefLevel bit=BitRefLevel(), Range *intersected=NULL)
MyVolumeFE feTRhs
calculate right hand side for tetrahedral elements
MyVolumeFE feMassRhs
calculate right hand side for tetrahedral elements
MyVolumeFE feTLhs
calculate left hand side for tetrahedral elements
MoFEMErrorCode setConvectiveMassOperators(string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool ale=false, bool linear=false)
boost::ptr_vector< MethodForForceScaling > methodsOp
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
MyVolumeFE feVelLhs
calculate left hand side for tetrahedral elements
static MoFEMErrorCode applyScale(const FEMethod *fe, boost::ptr_vector< MethodForForceScaling > &methods_op, VectorDouble &nf)
Managing BitRefLevels.
Body force data structure.
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode problem_basic_method_preProcess(const Problem *problem_ptr, BasicMethod &method, int verb=DEFAULT_VERBOSITY)=0
Set data for BasicMethod.
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)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
get dofs values
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1FieldData()
Return FTensor of rank 1, i.e. vector from filed data coeffinects.
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
const VectorInt & getIndices() const
Get global indices of dofs on entity.
structure to get information form mofem into EntitiesFieldData
Elastic material data structure.
intrusive_ptr for managing petsc objects
data structure for TS (time stepping) context
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Vector manager is used to create vectors \mofem_vectors.
Definition: VecManager.hpp:23