v0.15.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>
13#include <MethodForForceScaling.hpp>
14#include <DirichletBC.hpp>
15#include <MethodForForceScaling.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_NULLPTR), F(PETSC_NULLPTR) {
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>();
119 FTensor::Index<'i', 3> i;
120 FTensor::Index<'j', 3> j;
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 SETERRQ(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
389
391ConvectiveMassElement::OpMassRhs::doWork(int row_side, EntityType row_type,
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
411 FTensor::Index<'i', 3> i;
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)
459 FTensor::Index<'i', 3> i;
460 FTensor::Index<'j', 3> j;
461 FTensor::Index<'k', 3> k;
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, MOFEM_DATA_INCONSISTENCY, "%s",
556 ss.str().c_str());
557 }
558
559 FTensor::Index<'i', 3> i;
560 FTensor::Index<'j', 3> j;
561
562 {
563 int dd1 = 0;
564 // integrate element stiffness matrix
565 for (; dd1 < nb_row / 3; dd1++) {
567 &jac(0, 0), &jac(0, 1), &jac(0, 2), &jac(1, 0), &jac(1, 1),
568 &jac(1, 2), &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
569 for (int dd2 = 0; dd2 < nb_col / 3; dd2++) {
571 &k(3 * dd1 + 0, 3 * dd2 + 0), &k(3 * dd1 + 0, 3 * dd2 + 1),
572 &k(3 * dd1 + 0, 3 * dd2 + 2), &k(3 * dd1 + 1, 3 * dd2 + 0),
573 &k(3 * dd1 + 1, 3 * dd2 + 1), &k(3 * dd1 + 1, 3 * dd2 + 2),
574 &k(3 * dd1 + 2, 3 * dd2 + 0), &k(3 * dd1 + 2, 3 * dd2 + 1),
575 &k(3 * dd1 + 2, 3 * dd2 + 2));
576 t_k(i, j) += base * t_jac(i, j);
577 ++t_jac;
578 }
579 ++base;
580 // for(int rr1 = 0;rr1<3;rr1++) {
581 // for(int dd2 = 0;dd2<nb_col;dd2++) {
582 // k(3*dd1+rr1,dd2) += row_data.getN()(gg,dd1)*jac(rr1,dd2);
583 // }
584 // }
585 }
586 for (; dd1 != nb_base_functions; dd1++) {
587 ++base;
588 }
589 }
590 }
591
592 if (!forcesOnlyOnEntities.empty()) {
593 VectorInt indices = row_data.getIndices();
594 VectorDofs &dofs = row_data.getFieldDofs();
595 VectorDofs::iterator dit = dofs.begin();
596 for (int ii = 0; dit != dofs.end(); dit++, ii++) {
597 if (forcesOnlyOnEntities.find((*dit)->getEnt()) ==
598 forcesOnlyOnEntities.end()) {
599 indices[ii] = -1;
600 }
601 }
602 CHKERR MatSetValues(getFEMethod()->ts_B, nb_row, &indices[0], nb_col,
603 &col_data.getIndices()[0], &k(0, 0), ADD_VALUES);
604 } else {
605 CHKERR MatSetValues(getFEMethod()->ts_B, nb_row,
606 &row_data.getIndices()[0], nb_col,
607 &col_data.getIndices()[0], &k(0, 0), ADD_VALUES);
608 }
609 }
611}
612
614 const std::string field_name, const std::string col_field, BlockData &data,
615 CommonData &common_data)
616 : OpMassLhs_dM_dv(field_name, col_field, data, common_data) {}
617
619 EntitiesFieldData::EntData &col_data, int gg) {
621 FTensor::Index<'i', 3> i;
622 FTensor::Index<'j', 3> j;
623 FTensor::Index<'k', 3> k;
624 int nb_col = col_data.getIndices().size();
625 jac.clear();
626 FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
627 &jac(1, 0), &jac(1, 1), &jac(1, 2),
628 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
629 const int s = 3;
631 // T* d000, T* d001, T* d002,
632 // T* d010, T* d011, T* d012,
633 // T* d020, T* d021, T* d022,
634 // T* d100, T* d101, T* d102,
635 // T* d110, T* d111, T* d112,
636 // T* d120, T* d121, T* d122,
637 // T* d200, T* d201, T* d202,
638 // T* d210, T* d211, T* d212,
639 // T* d220, T* d221, T* d222,
640 &commonData.jacMass[gg](0, s + 0), &commonData.jacMass[gg](0, s + 1),
641 &commonData.jacMass[gg](0, s + 2), &commonData.jacMass[gg](0, s + 3),
642 &commonData.jacMass[gg](0, s + 4), &commonData.jacMass[gg](0, s + 5),
643 &commonData.jacMass[gg](0, s + 6), &commonData.jacMass[gg](0, s + 7),
644 &commonData.jacMass[gg](0, s + 8), &commonData.jacMass[gg](1, s + 0),
645 &commonData.jacMass[gg](1, s + 1), &commonData.jacMass[gg](1, s + 2),
646 &commonData.jacMass[gg](1, s + 3), &commonData.jacMass[gg](1, s + 4),
647 &commonData.jacMass[gg](1, s + 5), &commonData.jacMass[gg](1, s + 6),
648 &commonData.jacMass[gg](1, s + 7), &commonData.jacMass[gg](1, s + 8),
649 &commonData.jacMass[gg](2, s + 0), &commonData.jacMass[gg](2, s + 1),
650 &commonData.jacMass[gg](2, s + 2), &commonData.jacMass[gg](2, s + 3),
651 &commonData.jacMass[gg](2, s + 4), &commonData.jacMass[gg](2, s + 5),
652 &commonData.jacMass[gg](2, s + 6), &commonData.jacMass[gg](2, s + 7),
653 &commonData.jacMass[gg](2, s + 8));
654 double *diff_ptr =
655 const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
656 FTensor::Tensor1<double *, 3> diff(diff_ptr, &diff_ptr[1], &diff_ptr[2], 3);
657 for (int dd = 0; dd < nb_col / 3; dd++) {
658 t_jac(i, j) += t_mass3(i, j, k) * diff(k);
659 ++diff;
660 ++t_jac;
661 }
663}
664
666 const std::string field_name, const std::string col_field, BlockData &data,
667 CommonData &common_data)
668 : OpMassLhs_dM_dv(field_name, col_field, data, common_data) {}
669
671 EntitiesFieldData::EntData &col_data, int gg) {
673 int nb_col = col_data.getIndices().size();
674 jac.clear();
675 double *base_ptr = const_cast<double *>(&col_data.getN(gg)[0]);
676 FTensor::Tensor0<double *> base(base_ptr, 1);
677 double *diff_ptr =
678 const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
679 FTensor::Tensor1<double *, 3> diff(diff_ptr, &diff_ptr[1], &diff_ptr[2], 3);
680 FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
681 &jac(1, 0), &jac(1, 1), &jac(1, 2),
682 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
683 const int u = 3 + 9 + 9;
685 &commonData.jacMass[gg](0, u + 0), &commonData.jacMass[gg](0, u + 1),
686 &commonData.jacMass[gg](0, u + 2), &commonData.jacMass[gg](1, u + 0),
687 &commonData.jacMass[gg](1, u + 1), &commonData.jacMass[gg](1, u + 2),
688 &commonData.jacMass[gg](2, u + 0), &commonData.jacMass[gg](2, u + 1),
689 &commonData.jacMass[gg](2, u + 2));
690 const int s = 3 + 9 + 9 + 3;
692 // T* d000, T* d001, T* d002,
693 // T* d010, T* d011, T* d012,
694 // T* d020, T* d021, T* d022,
695 // T* d100, T* d101, T* d102,
696 // T* d110, T* d111, T* d112,
697 // T* d120, T* d121, T* d122,
698 // T* d200, T* d201, T* d202,
699 // T* d210, T* d211, T* d212,
700 // T* d220, T* d221, T* d222,
701 &commonData.jacMass[gg](0, s + 0), &commonData.jacMass[gg](0, s + 1),
702 &commonData.jacMass[gg](0, s + 2), &commonData.jacMass[gg](0, s + 3),
703 &commonData.jacMass[gg](0, s + 4), &commonData.jacMass[gg](0, s + 5),
704 &commonData.jacMass[gg](0, s + 6), &commonData.jacMass[gg](0, s + 7),
705 &commonData.jacMass[gg](0, s + 8), &commonData.jacMass[gg](1, s + 0),
706 &commonData.jacMass[gg](1, s + 1), &commonData.jacMass[gg](1, s + 2),
707 &commonData.jacMass[gg](1, s + 3), &commonData.jacMass[gg](1, s + 4),
708 &commonData.jacMass[gg](1, s + 5), &commonData.jacMass[gg](1, s + 6),
709 &commonData.jacMass[gg](1, s + 7), &commonData.jacMass[gg](1, s + 8),
710 &commonData.jacMass[gg](2, s + 0), &commonData.jacMass[gg](2, s + 1),
711 &commonData.jacMass[gg](2, s + 2), &commonData.jacMass[gg](2, s + 3),
712 &commonData.jacMass[gg](2, s + 4), &commonData.jacMass[gg](2, s + 5),
713 &commonData.jacMass[gg](2, s + 6), &commonData.jacMass[gg](2, s + 7),
714 &commonData.jacMass[gg](2, s + 8));
715 FTensor::Index<'i', 3> i;
716 FTensor::Index<'j', 3> j;
717 FTensor::Index<'k', 3> k;
718 for (int dd = 0; dd < nb_col / 3; dd++) {
719 t_jac(i, j) += t_mass1(i, j) * base * getFEMethod()->ts_a;
720 t_jac(i, j) += t_mass3(i, j, k) * diff(k);
721 ++base_ptr;
722 ++diff_ptr;
723 ++t_jac;
724 }
726}
727
729 BlockData &data,
730 CommonData &common_data,
734 dAta(data), commonData(common_data), V(v, true),
735 lInear(commonData.lInear) {}
736
738ConvectiveMassElement::OpEnergy::doWork(int row_side, EntityType row_type,
739 EntitiesFieldData::EntData &row_data) {
741
742 if (row_type != MBVERTEX) {
744 }
745 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
746 dAta.tEts.end()) {
748 }
749
750 {
751 double energy = 0;
752 for (unsigned int gg = 0; gg < row_data.getN().size1(); gg++) {
753 double val = getVolume() * getGaussPts()(3, gg);
754 double rho0 = dAta.rho0;
755 double rho;
756 if (lInear) {
757 rho = rho0;
758 } else {
759 h.resize(3, 3);
760 noalias(h) =
763 .size() > 0) {
764 H.resize(3, 3);
765 noalias(H) =
767 auto detH = determinantTensor3by3(H);
768 invH.resize(3, 3);
769 CHKERR invertTensor3by3(H, detH, invH);
770 F.resize(3, 3);
771 noalias(F) = prod(h, invH);
772 } else {
773 F.resize(3, 3);
774 noalias(F) = h;
775 }
776 double detF = determinantTensor3by3(F);
777 rho = detF * rho0;
778 }
779 v.resize(3);
781 energy += 0.5 * (rho * val) * inner_prod(v, v);
782 }
783 CHKERR VecSetValue(V, 0, energy, ADD_VALUES);
784 }
785
787}
788
790 const std::string field_name, BlockData &data, CommonData &common_data,
791 int tag, bool jacobian)
794 dAta(data), commonData(common_data), tAg(tag), jAcobian(jacobian),
795 fieldDisp(false) {}
796
798 int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data) {
800
801 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
802 dAta.tEts.end()) {
804 }
805
806 // do it only once, no need to repeat this for edges,faces or tets
807 if (row_type != MBVERTEX)
809
810 int nb_dofs = row_data.getIndices().size();
811 if (nb_dofs == 0)
813
814 {
815
816 v.resize(3);
817 dot_w.resize(3);
818 h.resize(3, 3);
819 h.clear();
820 F.resize(3, 3);
821 dot_W.resize(3);
822 dot_W.clear();
823 H.resize(3, 3);
824 H.clear();
825 invH.resize(3, 3);
826 invH.clear();
827 dot_u.resize(3);
828 for (int dd = 0; dd < 3; dd++) {
829 H(dd, dd) = 1;
830 invH(dd, dd) = 1;
831 }
832
833 a_res.resize(3);
834 int nb_gauss_pts = row_data.getN().size1();
835 commonData.valVel.resize(nb_gauss_pts);
836 commonData.jacVelRowPtr.resize(nb_gauss_pts);
837 commonData.jacVel.resize(nb_gauss_pts);
838
839 int nb_active_vars = 0;
840 for (int gg = 0; gg < nb_gauss_pts; gg++) {
841
842 if (gg == 0) {
843
844 trace_on(tAg);
845
846 for (int nn1 = 0; nn1 < 3; nn1++) { // 0
847 v[nn1] <<=
849 nb_active_vars++;
850 }
851 for (int nn1 = 0; nn1 < 3; nn1++) { // 3
852 dot_w[nn1] <<=
854 [gg][nn1];
855 nb_active_vars++;
856 }
858 .size() > 0) {
859 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3 = 6
860 for (int nn2 = 0; nn2 < 3; nn2++) {
861 h(nn1, nn2) <<=
863 nn1, nn2);
864 if (fieldDisp) {
865 if (nn1 == nn2) {
866 h(nn1, nn2) += 1;
867 }
868 }
869 nb_active_vars++;
870 }
871 }
872 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9
873 dot_W[nn1] <<=
875 .dataAtGaussPts["DOT_" + commonData.meshPositions][gg][nn1];
876 nb_active_vars++;
877 }
878 }
880 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9+3
881 for (int nn2 = 0; nn2 < 3; nn2++) {
882 H(nn1, nn2) <<=
884 nn2);
885 nb_active_vars++;
886 }
887 }
888 }
889 detH = 1;
890
891 FTensor::Index<'i', 3> i;
892 FTensor::Index<'j', 3> j;
893 FTensor::Index<'k', 3> k;
894
898 auto t_invH = getFTensor2FromArray3by3(invH, FTensor::Number<0>(), 0);
899 auto t_dot_u =
900 FTensor::Tensor1<adouble *, 3>{&dot_u[0], &dot_u[1], &dot_u[2]};
901 auto t_dot_w =
902 FTensor::Tensor1<adouble *, 3>{&dot_w[0], &dot_w[1], &dot_w[2]};
903 auto t_dot_W =
904 FTensor::Tensor1<adouble *, 3>{&dot_W[0], &dot_W[1], &dot_W[2]};
905 auto t_v = FTensor::Tensor1<adouble *, 3>{&v[0], &v[1], &v[2]};
906 auto t_a_res =
907 FTensor::Tensor1<adouble *, 3>{&a_res[0], &a_res[1], &a_res[2]};
908
910 detH = determinantTensor3by3(H);
911 CHKERR invertTensor3by3(H, detH, invH);
912 t_F(i, j) = t_h(i, k) * t_invH(k, j);
913 } else {
914 t_F(i, j) = t_h(i, j);
915 }
916
917 t_dot_u(i) = t_dot_w(i) + t_F(i, j) * t_dot_W(j);
918 t_a_res(i) = t_v(i) - t_dot_u(i);
919 t_a_res(i) *= detH;
920
921 // dependant
922 VectorDouble &res = commonData.valVel[gg];
923 res.resize(3);
924 for (int rr = 0; rr < 3; rr++) {
925 a_res[rr] >>= res[rr];
926 }
927 trace_off();
928 }
929
930 active.resize(nb_active_vars);
931 int aa = 0;
932 for (int nn1 = 0; nn1 < 3; nn1++) {
933 active[aa++] =
935 }
936 for (int nn1 = 0; nn1 < 3; nn1++) {
937 active[aa++] =
939 .dataAtGaussPts["DOT_" + commonData.spatialPositions][gg][nn1];
940 }
942 0) {
943 for (int nn1 = 0; nn1 < 3; nn1++) {
944 for (int nn2 = 0; nn2 < 3; nn2++) {
945 if (fieldDisp && nn1 == nn2) {
946 active[aa++] =
948 nn1, nn2) +
949 1;
950 } else {
951 active[aa++] =
953 nn1, nn2);
954 }
955 }
956 }
957 for (int nn1 = 0; nn1 < 3; nn1++) {
958 active[aa++] =
960 .dataAtGaussPts["DOT_" + commonData.meshPositions][gg][nn1];
961 }
962 }
964 for (int nn1 = 0; nn1 < 3; nn1++) {
965 for (int nn2 = 0; nn2 < 3; nn2++) {
966 active[aa++] =
968 nn2);
969 }
970 }
971 }
972
973 if (!jAcobian) {
974 VectorDouble &res = commonData.valVel[gg];
975 if (gg > 0) {
976 res.resize(3);
977 int r;
978 r = ::function(tAg, 3, nb_active_vars, &active[0], &res[0]);
979 if (r != 3) {
980 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
981 "ADOL-C function evaluation with error");
982 }
983 }
984 double val = getVolume() * getGaussPts()(3, gg);
985 res *= val;
986 } else {
987 commonData.jacVelRowPtr[gg].resize(3);
988 commonData.jacVel[gg].resize(3, nb_active_vars);
989 for (int nn1 = 0; nn1 < 3; nn1++) {
990 (commonData.jacVelRowPtr[gg])[nn1] = &(commonData.jacVel[gg](nn1, 0));
991 }
992 int r;
993 r = jacobian(tAg, 3, nb_active_vars, &active[0],
994 &(commonData.jacVelRowPtr[gg])[0]);
995 if (r != 3) {
996 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
997 "ADOL-C function evaluation with error");
998 }
999 double val = getVolume() * getGaussPts()(3, gg);
1000 commonData.jacVel[gg] *= val;
1001 // std::cerr << gg << " : " << commonData.jacVel[gg] << std::endl;
1002 }
1003 }
1004 }
1006}
1007
1013
1015 int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data) {
1017
1018 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
1019 dAta.tEts.end()) {
1021 }
1022 int nb_dofs = row_data.getIndices().size();
1023 if (nb_dofs == 0)
1025
1026 auto base = row_data.getFTensor0N();
1027 int nb_base_functions = row_data.getN().size2();
1028 FTensor::Index<'i', 3> i;
1029
1030 {
1031
1032 nf.resize(nb_dofs);
1033 nf.clear();
1034
1035 for (unsigned int gg = 0; gg < row_data.getN().size1(); gg++) {
1036 FTensor::Tensor1<double *, 3> t_nf(&nf[0], &nf[1], &nf[2], 3);
1038 &commonData.valVel[gg][1],
1039 &commonData.valVel[gg][2]);
1040 int dd = 0;
1041 for (; dd < nb_dofs / 3; dd++) {
1042 t_nf(i) += base * res(i);
1043 ++base;
1044 ++t_nf;
1045 }
1046 for (; dd != nb_base_functions; dd++) {
1047 ++base;
1048 }
1049 }
1050
1051 if (row_data.getIndices().size() > 3 * row_data.getN().size2()) {
1052 SETERRQ(PETSC_COMM_SELF, 1, "data inconsistency");
1053 }
1054 CHKERR VecSetValues(getFEMethod()->ts_F, row_data.getIndices().size(),
1055 &row_data.getIndices()[0], &nf[0], ADD_VALUES);
1056 }
1058}
1059
1061 const std::string vel_field, const std::string field_name, BlockData &data,
1062 CommonData &common_data)
1063 : OpMassLhs_dM_dv(vel_field, field_name, data, common_data) {}
1064
1066 EntitiesFieldData::EntData &col_data, int gg) {
1068 int nb_col = col_data.getIndices().size();
1069 jac.clear();
1070 if (!nb_col)
1072 double *base_ptr = const_cast<double *>(&col_data.getN(gg)[0]);
1073 FTensor::Tensor0<double *> base(base_ptr, 1);
1074 FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
1075 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1076 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1078 &commonData.jacVel[gg](0, 0), &commonData.jacVel[gg](0, 1),
1079 &commonData.jacVel[gg](0, 2), &commonData.jacVel[gg](1, 0),
1080 &commonData.jacVel[gg](1, 1), &commonData.jacVel[gg](1, 2),
1081 &commonData.jacVel[gg](2, 0), &commonData.jacVel[gg](2, 1),
1082 &commonData.jacVel[gg](2, 2));
1083 FTensor::Index<'i', 3> i;
1084 FTensor::Index<'j', 3> j;
1085 for (int dd = 0; dd < nb_col / 3; dd++) {
1086 t_jac(i, j) += t_mass1(i, j) * base;
1087 ++base_ptr;
1088 ++t_jac;
1089 }
1090
1092}
1093
1095 const std::string vel_field, const std::string field_name, BlockData &data,
1096 CommonData &common_data)
1097 : OpVelocityLhs_dV_dv(vel_field, field_name, data, common_data) {}
1098
1100 EntitiesFieldData::EntData &col_data, int gg) {
1102 int nb_col = col_data.getIndices().size();
1103 jac.clear();
1104 if (!nb_col)
1106 double *base_ptr = const_cast<double *>(&col_data.getN(gg)[0]);
1107 FTensor::Tensor0<double *> base(base_ptr, 1);
1108 FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
1109 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1110 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1111 const int u = 3;
1113 &commonData.jacVel[gg](0, u + 0), &commonData.jacVel[gg](0, u + 1),
1114 &commonData.jacVel[gg](0, u + 2), &commonData.jacVel[gg](1, u + 0),
1115 &commonData.jacVel[gg](1, u + 1), &commonData.jacVel[gg](1, u + 2),
1116 &commonData.jacVel[gg](2, u + 0), &commonData.jacVel[gg](2, u + 1),
1117 &commonData.jacVel[gg](2, u + 2));
1118 FTensor::Index<'i', 3> i;
1119 FTensor::Index<'j', 3> j;
1120 FTensor::Index<'k', 3> k;
1121 if (commonData.dataAtGaussPts["DOT_" + commonData.meshPositions].size() ==
1122 0) {
1123
1124 for (int dd = 0; dd < nb_col / 3; dd++) {
1125 t_jac(i, j) += t_mass1(i, j) * base * getFEMethod()->ts_a;
1126 ++base_ptr;
1127 ++t_jac;
1128 }
1129 } else {
1130 double *diff_ptr =
1131 const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
1132 FTensor::Tensor1<double *, 3> diff(diff_ptr, &diff_ptr[1], &diff_ptr[2], 3);
1133 const int s = 3 + 3;
1135 // T* d000, T* d001, T* d002,
1136 // T* d010, T* d011, T* d012,
1137 // T* d020, T* d021, T* d022,
1138 // T* d100, T* d101, T* d102,
1139 // T* d110, T* d111, T* d112,
1140 // T* d120, T* d121, T* d122,
1141 // T* d200, T* d201, T* d202,
1142 // T* d210, T* d211, T* d212,
1143 // T* d220, T* d221, T* d222,
1144 &commonData.jacVel[gg](0, s + 0), &commonData.jacVel[gg](0, s + 1),
1145 &commonData.jacVel[gg](0, s + 2), &commonData.jacVel[gg](0, s + 3),
1146 &commonData.jacVel[gg](0, s + 4), &commonData.jacVel[gg](0, s + 5),
1147 &commonData.jacVel[gg](0, s + 6), &commonData.jacVel[gg](0, s + 7),
1148 &commonData.jacVel[gg](0, s + 8), &commonData.jacVel[gg](1, s + 0),
1149 &commonData.jacVel[gg](1, s + 1), &commonData.jacVel[gg](1, s + 2),
1150 &commonData.jacVel[gg](1, s + 3), &commonData.jacVel[gg](1, s + 4),
1151 &commonData.jacVel[gg](1, s + 5), &commonData.jacVel[gg](1, s + 6),
1152 &commonData.jacVel[gg](1, s + 7), &commonData.jacVel[gg](1, s + 8),
1153 &commonData.jacVel[gg](2, s + 0), &commonData.jacVel[gg](2, s + 1),
1154 &commonData.jacVel[gg](2, s + 2), &commonData.jacVel[gg](2, s + 3),
1155 &commonData.jacVel[gg](2, s + 4), &commonData.jacVel[gg](2, s + 5),
1156 &commonData.jacVel[gg](2, s + 6), &commonData.jacVel[gg](2, s + 7),
1157 &commonData.jacVel[gg](2, s + 8));
1158 for (int dd = 0; dd < nb_col / 3; dd++) {
1159 t_jac(i, j) += t_mass1(i, j) * base * getFEMethod()->ts_a;
1160 t_jac(i, j) += t_mass3(i, j, k) * diff(k);
1161 ++base_ptr;
1162 ++diff_ptr;
1163 ++t_jac;
1164 }
1165 }
1167}
1168
1170 const std::string vel_field, const std::string field_name, BlockData &data,
1171 CommonData &common_data)
1172 : OpVelocityLhs_dV_dv(vel_field, field_name, data, common_data) {}
1173
1175 EntitiesFieldData::EntData &col_data, int gg) {
1177 int nb_col = col_data.getIndices().size();
1178 jac.clear();
1179 if (!nb_col)
1181 double *base_ptr = const_cast<double *>(&col_data.getN(gg)[0]);
1182 FTensor::Tensor0<double *> base(base_ptr, 1);
1183 double *diff_ptr =
1184 const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
1185 FTensor::Tensor1<double *, 3> diff(diff_ptr, &diff_ptr[1], &diff_ptr[2], 3);
1186 FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
1187 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1188 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1189 const int u = 3 + 3 + 9;
1191 &commonData.jacVel[gg](0, u + 0), &commonData.jacVel[gg](0, u + 1),
1192 &commonData.jacVel[gg](0, u + 2), &commonData.jacVel[gg](1, u + 0),
1193 &commonData.jacVel[gg](1, u + 1), &commonData.jacVel[gg](1, u + 2),
1194 &commonData.jacVel[gg](2, u + 0), &commonData.jacVel[gg](2, u + 1),
1195 &commonData.jacVel[gg](2, u + 2));
1196 const int s = 3 + 3 + 9 + 3;
1198 // T* d000, T* d001, T* d002,
1199 // T* d010, T* d011, T* d012,
1200 // T* d020, T* d021, T* d022,
1201 // T* d100, T* d101, T* d102,
1202 // T* d110, T* d111, T* d112,
1203 // T* d120, T* d121, T* d122,
1204 // T* d200, T* d201, T* d202,
1205 // T* d210, T* d211, T* d212,
1206 // T* d220, T* d221, T* d222,
1207 &commonData.jacVel[gg](0, s + 0), &commonData.jacVel[gg](0, s + 1),
1208 &commonData.jacVel[gg](0, s + 2), &commonData.jacVel[gg](0, s + 3),
1209 &commonData.jacVel[gg](0, s + 4), &commonData.jacVel[gg](0, s + 5),
1210 &commonData.jacVel[gg](0, s + 6), &commonData.jacVel[gg](0, s + 7),
1211 &commonData.jacVel[gg](0, s + 8), &commonData.jacVel[gg](1, s + 0),
1212 &commonData.jacVel[gg](1, s + 1), &commonData.jacVel[gg](1, s + 2),
1213 &commonData.jacVel[gg](1, s + 3), &commonData.jacVel[gg](1, s + 4),
1214 &commonData.jacVel[gg](1, s + 5), &commonData.jacVel[gg](1, s + 6),
1215 &commonData.jacVel[gg](1, s + 7), &commonData.jacVel[gg](1, s + 8),
1216 &commonData.jacVel[gg](2, s + 0), &commonData.jacVel[gg](2, s + 1),
1217 &commonData.jacVel[gg](2, s + 2), &commonData.jacVel[gg](2, s + 3),
1218 &commonData.jacVel[gg](2, s + 4), &commonData.jacVel[gg](2, s + 5),
1219 &commonData.jacVel[gg](2, s + 6), &commonData.jacVel[gg](2, s + 7),
1220 &commonData.jacVel[gg](2, s + 8));
1221 FTensor::Index<'i', 3> i;
1222 FTensor::Index<'j', 3> j;
1223 FTensor::Index<'k', 3> k;
1224 for (int dd = 0; dd < nb_col / 3; dd++) {
1225 t_jac(i, j) += t_mass1(i, j) * base * getFEMethod()->ts_a;
1226 t_jac(i, j) += t_mass3(i, j, k) * diff(k);
1227 ++base_ptr;
1228 ++diff_ptr;
1229 ++t_jac;
1230 }
1232}
1233
1236 BlockData &data,
1237 CommonData &common_data, int tag,
1238 bool jacobian)
1241 dAta(data), commonData(common_data), tAg(tag), jAcobian(jacobian),
1242 fieldDisp(false) {}
1243
1246 int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data) {
1248
1249 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
1250 dAta.tEts.end()) {
1252 }
1253
1254 // do it only once, no need to repeat this for edges,faces or tets
1255 if (row_type != MBVERTEX)
1257
1258 int nb_dofs = row_data.getIndices().size();
1259 if (nb_dofs == 0)
1261
1262 try {
1263
1264 a.resize(3);
1265 v.resize(3);
1266 g.resize(3, 3);
1267 G.resize(3, 3);
1268 h.resize(3, 3);
1269 F.resize(3, 3);
1270 H.resize(3, 3);
1271 H.clear();
1272 invH.resize(3, 3);
1273 invH.clear();
1274 for (int dd = 0; dd < 3; dd++) {
1275 H(dd, dd) = 1;
1276 invH(dd, dd) = 1;
1277 }
1278
1279 int nb_gauss_pts = row_data.getN().size1();
1280 commonData.valT.resize(nb_gauss_pts);
1281 commonData.jacTRowPtr.resize(nb_gauss_pts);
1282 commonData.jacT.resize(nb_gauss_pts);
1283
1284 int nb_active_vars = 0;
1285 for (int gg = 0; gg < nb_gauss_pts; gg++) {
1286
1287 if (gg == 0) {
1288
1289 trace_on(tAg);
1290
1291 for (int nn1 = 0; nn1 < 3; nn1++) { // 0
1292 a[nn1] <<=
1294 [gg][nn1];
1295 nb_active_vars++;
1296 }
1297
1298 for (int nn1 = 0; nn1 < 3; nn1++) { // 3
1299 v[nn1] <<=
1301 nb_active_vars++;
1302 }
1303 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3
1304 for (int nn2 = 0; nn2 < 3; nn2++) {
1305 g(nn1, nn2) <<=
1307 nn1, nn2);
1308 nb_active_vars++;
1309 }
1310 }
1311 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9
1312 for (int nn2 = 0; nn2 < 3; nn2++) {
1313 h(nn1, nn2) <<=
1315 nn2);
1316 nb_active_vars++;
1317 if (fieldDisp) {
1318 if (nn1 == nn2) {
1319 h(nn1, nn2) += 1;
1320 }
1321 }
1322 }
1323 }
1325 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9+9
1326 for (int nn2 = 0; nn2 < 3; nn2++) {
1327 H(nn1, nn2) <<=
1329 nn2);
1330 nb_active_vars++;
1331 }
1332 }
1333 }
1334 adouble detH;
1335 detH = 1;
1337 detH = determinantTensor3by3(H);
1338 }
1339 CHKERR invertTensor3by3(H, detH, invH);
1340
1341 FTensor::Index<'i', 3> i;
1342 FTensor::Index<'j', 3> j;
1343 FTensor::Index<'k', 3> k;
1344
1345 a_T.resize(3);
1346
1348 auto t_invH = getFTensor2FromArray3by3(invH, FTensor::Number<0>(), 0);
1351 auto t_G = getFTensor2FromArray3by3(G, FTensor::Number<0>(), 0);
1352
1353 auto t_a = FTensor::Tensor1<adouble *, 3>{&a[0], &a[1], &a[2]};
1354 auto t_v = FTensor::Tensor1<adouble *, 3>{&v[0], &v[1], &v[2]};
1355 auto t_a_T = FTensor::Tensor1<adouble *, 3>{&a_T[0], &a_T[1], &a_T[2]};
1356
1357 t_F(i, j) = t_h(i, k) * t_invH(k, j);
1358 t_G(i, j) = t_g(i, k) * t_invH(k, j);
1359 t_a_T(i) = t_F(k, i) * t_a(k) + t_G(k, i) * t_v(k);
1360 const auto rho0 = dAta.rho0;
1361 t_a_T(i) = -rho0 * detH;
1362
1363 commonData.valT[gg].resize(3);
1364 for (int nn = 0; nn < 3; nn++) {
1365 a_T[nn] >>= (commonData.valT[gg])[nn];
1366 }
1367 trace_off();
1368 }
1369
1370 active.resize(nb_active_vars);
1371 int aa = 0;
1372 for (int nn1 = 0; nn1 < 3; nn1++) { // 0
1373 active[aa++] =
1375 .dataAtGaussPts["DOT_" + commonData.spatialVelocities][gg][nn1];
1376 }
1377
1378 for (int nn1 = 0; nn1 < 3; nn1++) { // 3
1379 active[aa++] =
1381 }
1382 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3
1383 for (int nn2 = 0; nn2 < 3; nn2++) {
1384 active[aa++] =
1386 nn2);
1387 }
1388 }
1389 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9
1390 for (int nn2 = 0; nn2 < 3; nn2++) {
1391 if (fieldDisp && nn1 == nn2) {
1392 active[aa++] =
1394 nn1, nn2) +
1395 1;
1396 } else {
1397 active[aa++] =
1399 nn2);
1400 }
1401 }
1402 }
1404 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9+9
1405 for (int nn2 = 0; nn2 < 3; nn2++) {
1406 active[aa++] =
1408 nn2);
1409 }
1410 }
1411 }
1412
1413 if (!jAcobian) {
1414 VectorDouble &res = commonData.valT[gg];
1415 if (gg > 0) {
1416 res.resize(3);
1417 int r;
1418 r = ::function(tAg, 3, nb_active_vars, &active[0], &res[0]);
1419 if (r != 3) { // function is locally analytic
1420 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1421 "ADOL-C function evaluation with error r = %d", r);
1422 }
1423 }
1424 double val = getVolume() * getGaussPts()(3, gg);
1425 res *= val;
1426 } else {
1427 commonData.jacTRowPtr[gg].resize(3);
1428 commonData.jacT[gg].resize(3, nb_active_vars);
1429 for (int nn1 = 0; nn1 < 3; nn1++) {
1430 (commonData.jacTRowPtr[gg])[nn1] = &(commonData.jacT[gg](nn1, 0));
1431 }
1432 int r;
1433 r = jacobian(tAg, 3, nb_active_vars, &active[0],
1434 &(commonData.jacTRowPtr[gg])[0]);
1435 if (r != 3) {
1436 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1437 "ADOL-C function evaluation with error");
1438 }
1439 double val = getVolume() * getGaussPts()(3, gg);
1440 commonData.jacT[gg] *= val;
1441 }
1442 }
1443
1444 } catch (const std::exception &ex) {
1445 std::ostringstream ss;
1446 ss << "throw in method: " << ex.what() << std::endl;
1447 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "%s", ss.str().c_str());
1448 }
1449
1451}
1452
1455 BlockData &data,
1456 CommonData &common_data,
1457 Range *forcesonlyonentities_ptr)
1460 dAta(data), commonData(common_data) {
1461 if (forcesonlyonentities_ptr != NULL) {
1462 forcesOnlyOnEntities = *forcesonlyonentities_ptr;
1463 }
1464}
1465
1468 int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data) {
1470
1471 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
1472 dAta.tEts.end()) {
1474 }
1475 int nb_dofs = row_data.getIndices().size();
1476 if (nb_dofs == 0)
1478
1479 try {
1480
1481 nf.resize(nb_dofs);
1482 nf.clear();
1483
1484 auto base = row_data.getFTensor0N();
1485 int nb_base_functions = row_data.getN().size2();
1486 FTensor::Index<'i', 3> i;
1487
1488 for (unsigned int gg = 0; gg < row_data.getN().size1(); gg++) {
1489 FTensor::Tensor1<double *, 3> t_nf(&nf[0], &nf[1], &nf[2], 3);
1491 &commonData.valT[gg][1],
1492 &commonData.valT[gg][2]);
1493 int dd = 0;
1494 for (; dd < nb_dofs / 3; dd++) {
1495 t_nf(i) += base * res(i);
1496 ++base;
1497 ++t_nf;
1498 }
1499 for (; dd != nb_base_functions; dd++) {
1500 ++base;
1501 }
1502 }
1503
1504 if (row_data.getIndices().size() > 3 * row_data.getN().size2()) {
1505 SETERRQ(PETSC_COMM_SELF, 1, "data inconsistency");
1506 }
1507 if (!forcesOnlyOnEntities.empty()) {
1508 VectorInt indices = row_data.getIndices();
1509 VectorDofs &dofs = row_data.getFieldDofs();
1510 VectorDofs::iterator dit = dofs.begin();
1511 for (int ii = 0; dit != dofs.end(); dit++, ii++) {
1512 if (forcesOnlyOnEntities.find((*dit)->getEnt()) ==
1513 forcesOnlyOnEntities.end()) {
1514 // std::cerr << **dit << std::endl;
1515 indices[ii] = -1;
1516 }
1517 }
1518 // std::cerr << indices << std::endl;
1519 CHKERR VecSetValues(getFEMethod()->ts_F, indices.size(), &indices[0],
1520 &nf[0], ADD_VALUES);
1521 } else {
1522 CHKERR VecSetValues(getFEMethod()->ts_F, row_data.getIndices().size(),
1523 &row_data.getIndices()[0], &nf[0], ADD_VALUES);
1524 }
1525
1526 } catch (const std::exception &ex) {
1527 std::ostringstream ss;
1528 ss << "throw in method: " << ex.what() << std::endl;
1529 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "%s", ss.str().c_str());
1530 }
1531
1533}
1534
1536 OpEshelbyDynamicMaterialMomentumLhs_dv(const std::string vel_field,
1537 const std::string field_name,
1538 BlockData &data,
1539 CommonData &common_data,
1540 Range *forcesonlyonentities_ptr)
1542 vel_field, field_name, data, common_data, forcesonlyonentities_ptr) {}
1543
1546 EntitiesFieldData::EntData &col_data, int gg) {
1548 int nb_col = col_data.getIndices().size();
1549 jac.clear();
1550 double *base_ptr = const_cast<double *>(&col_data.getN(gg)[0]);
1551 FTensor::Tensor0<double *> base(base_ptr, 1);
1552 double *diff_ptr =
1553 const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
1554 FTensor::Tensor1<double *, 3> diff(diff_ptr, &diff_ptr[1], &diff_ptr[2], 3);
1555 FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
1556 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1557 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1558 const int u = 3;
1560 &commonData.jacT[gg](0, u + 0), &commonData.jacT[gg](0, u + 1),
1561 &commonData.jacT[gg](0, u + 2), &commonData.jacT[gg](1, u + 0),
1562 &commonData.jacT[gg](1, u + 1), &commonData.jacT[gg](1, u + 2),
1563 &commonData.jacT[gg](2, u + 0), &commonData.jacT[gg](2, u + 1),
1564 &commonData.jacT[gg](2, u + 2));
1565 const int s = 3 + 3;
1567 // T* d000, T* d001, T* d002,
1568 // T* d010, T* d011, T* d012,
1569 // T* d020, T* d021, T* d022,
1570 // T* d100, T* d101, T* d102,
1571 // T* d110, T* d111, T* d112,
1572 // T* d120, T* d121, T* d122,
1573 // T* d200, T* d201, T* d202,
1574 // T* d210, T* d211, T* d212,
1575 // T* d220, T* d221, T* d222,
1576 &commonData.jacT[gg](0, s + 0), &commonData.jacT[gg](0, s + 1),
1577 &commonData.jacT[gg](0, s + 2), &commonData.jacT[gg](0, s + 3),
1578 &commonData.jacT[gg](0, s + 4), &commonData.jacT[gg](0, s + 5),
1579 &commonData.jacT[gg](0, s + 6), &commonData.jacT[gg](0, s + 7),
1580 &commonData.jacT[gg](0, s + 8), &commonData.jacT[gg](1, s + 0),
1581 &commonData.jacT[gg](1, s + 1), &commonData.jacT[gg](1, s + 2),
1582 &commonData.jacT[gg](1, s + 3), &commonData.jacT[gg](1, s + 4),
1583 &commonData.jacT[gg](1, s + 5), &commonData.jacT[gg](1, s + 6),
1584 &commonData.jacT[gg](1, s + 7), &commonData.jacT[gg](1, s + 8),
1585 &commonData.jacT[gg](2, s + 0), &commonData.jacT[gg](2, s + 1),
1586 &commonData.jacT[gg](2, s + 2), &commonData.jacT[gg](2, s + 3),
1587 &commonData.jacT[gg](2, s + 4), &commonData.jacT[gg](2, s + 5),
1588 &commonData.jacT[gg](2, s + 6), &commonData.jacT[gg](2, s + 7),
1589 &commonData.jacT[gg](2, s + 8));
1590 FTensor::Index<'i', 3> i;
1591 FTensor::Index<'j', 3> j;
1592 FTensor::Index<'k', 3> k;
1593 for (int dd = 0; dd < nb_col / 3; dd++) {
1594 t_jac(i, j) += t_mass1(i, j) * base * getFEMethod()->ts_a;
1595 t_jac(i, j) += t_mass3(i, j, k) * diff(k);
1596 ++base_ptr;
1597 ++diff_ptr;
1598 ++t_jac;
1599 }
1601}
1602
1604 OpEshelbyDynamicMaterialMomentumLhs_dx(const std::string vel_field,
1605 const std::string field_name,
1606 BlockData &data,
1607 CommonData &common_data,
1608 Range *forcesonlyonentities_ptr)
1610 vel_field, field_name, data, common_data, forcesonlyonentities_ptr) {}
1611
1614 EntitiesFieldData::EntData &col_data, int gg) {
1616 int nb_col = col_data.getIndices().size();
1617 jac.clear();
1618 double *diff_ptr =
1619 const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
1620 FTensor::Tensor1<double *, 3> diff(diff_ptr, &diff_ptr[1], &diff_ptr[2], 3);
1621 FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
1622 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1623 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1624 const int s = 3 + 3 + 9;
1626 // T* d000, T* d001, T* d002,
1627 // T* d010, T* d011, T* d012,
1628 // T* d020, T* d021, T* d022,
1629 // T* d100, T* d101, T* d102,
1630 // T* d110, T* d111, T* d112,
1631 // T* d120, T* d121, T* d122,
1632 // T* d200, T* d201, T* d202,
1633 // T* d210, T* d211, T* d212,
1634 // T* d220, T* d221, T* d222,
1635 &commonData.jacT[gg](0, s + 0), &commonData.jacT[gg](0, s + 1),
1636 &commonData.jacT[gg](0, s + 2), &commonData.jacT[gg](0, s + 3),
1637 &commonData.jacT[gg](0, s + 4), &commonData.jacT[gg](0, s + 5),
1638 &commonData.jacT[gg](0, s + 6), &commonData.jacT[gg](0, s + 7),
1639 &commonData.jacT[gg](0, s + 8), &commonData.jacT[gg](1, s + 0),
1640 &commonData.jacT[gg](1, s + 1), &commonData.jacT[gg](1, s + 2),
1641 &commonData.jacT[gg](1, s + 3), &commonData.jacT[gg](1, s + 4),
1642 &commonData.jacT[gg](1, s + 5), &commonData.jacT[gg](1, s + 6),
1643 &commonData.jacT[gg](1, s + 7), &commonData.jacT[gg](1, s + 8),
1644 &commonData.jacT[gg](2, s + 0), &commonData.jacT[gg](2, s + 1),
1645 &commonData.jacT[gg](2, s + 2), &commonData.jacT[gg](2, s + 3),
1646 &commonData.jacT[gg](2, s + 4), &commonData.jacT[gg](2, s + 5),
1647 &commonData.jacT[gg](2, s + 6), &commonData.jacT[gg](2, s + 7),
1648 &commonData.jacT[gg](2, s + 8));
1649 FTensor::Index<'i', 3> i;
1650 FTensor::Index<'j', 3> j;
1651 FTensor::Index<'k', 3> k;
1652 for (int dd = 0; dd < nb_col / 3; dd++) {
1653 t_jac(i, j) += t_mass3(i, j, k) * diff(k);
1654 ++diff_ptr;
1655 ++t_jac;
1656 }
1658}
1659
1661 OpEshelbyDynamicMaterialMomentumLhs_dX(const std::string vel_field,
1662 const std::string field_name,
1663 BlockData &data,
1664 CommonData &common_data,
1665 Range *forcesonlyonentities_ptr)
1667 vel_field, field_name, data, common_data, forcesonlyonentities_ptr) {}
1668
1671 EntitiesFieldData::EntData &col_data, int gg) {
1673 int nb_col = col_data.getIndices().size();
1674 jac.clear();
1675 double *diff_ptr =
1676 const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
1677 FTensor::Tensor1<double *, 3> diff(diff_ptr, &diff_ptr[1], &diff_ptr[2], 3);
1678 FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
1679 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1680 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1681 const int s = 3 + 3 + 9 + 9;
1683 // T* d000, T* d001, T* d002,
1684 // T* d010, T* d011, T* d012,
1685 // T* d020, T* d021, T* d022,
1686 // T* d100, T* d101, T* d102,
1687 // T* d110, T* d111, T* d112,
1688 // T* d120, T* d121, T* d122,
1689 // T* d200, T* d201, T* d202,
1690 // T* d210, T* d211, T* d212,
1691 // T* d220, T* d221, T* d222,
1692 &commonData.jacT[gg](0, s + 0), &commonData.jacT[gg](0, s + 1),
1693 &commonData.jacT[gg](0, s + 2), &commonData.jacT[gg](0, s + 3),
1694 &commonData.jacT[gg](0, s + 4), &commonData.jacT[gg](0, s + 5),
1695 &commonData.jacT[gg](0, s + 6), &commonData.jacT[gg](0, s + 7),
1696 &commonData.jacT[gg](0, s + 8), &commonData.jacT[gg](1, s + 0),
1697 &commonData.jacT[gg](1, s + 1), &commonData.jacT[gg](1, s + 2),
1698 &commonData.jacT[gg](1, s + 3), &commonData.jacT[gg](1, s + 4),
1699 &commonData.jacT[gg](1, s + 5), &commonData.jacT[gg](1, s + 6),
1700 &commonData.jacT[gg](1, s + 7), &commonData.jacT[gg](1, s + 8),
1701 &commonData.jacT[gg](2, s + 0), &commonData.jacT[gg](2, s + 1),
1702 &commonData.jacT[gg](2, s + 2), &commonData.jacT[gg](2, s + 3),
1703 &commonData.jacT[gg](2, s + 4), &commonData.jacT[gg](2, s + 5),
1704 &commonData.jacT[gg](2, s + 6), &commonData.jacT[gg](2, s + 7),
1705 &commonData.jacT[gg](2, s + 8));
1706 FTensor::Index<'i', 3> i;
1707 FTensor::Index<'j', 3> j;
1708 FTensor::Index<'k', 3> k;
1709 for (int dd = 0; dd < nb_col / 3; dd++) {
1710 t_jac(i, j) += t_mass3(i, j, k) * diff(k);
1711 ++diff_ptr;
1712 ++t_jac;
1713 }
1715}
1716
1718 MoFEM::Interface &m_field, TS _ts, const std::string velocity_field,
1719 const std::string spatial_position_field)
1720 : mField(m_field), tS(_ts), velocityField(velocity_field),
1721 spatialPositionField(spatial_position_field), jacobianLag(-1) {}
1722
1725
1726 switch (ts_ctx) {
1727 case CTX_TSSETIFUNCTION: {
1728 snes_f = ts_F;
1729 // FIXME: This global scattering because Kuu problem and Dynamic problem
1730 // not share partitions. Both problem should use the same partitioning to
1731 // resolve this problem.
1732 CHKERR mField.getInterface<VecManager>()->setGlobalGhostVector(
1733 problemPtr, COL, ts_u, INSERT_VALUES, SCATTER_REVERSE);
1734 CHKERR mField.getInterface<VecManager>()->setOtherGlobalGhostVector(
1735 problemPtr, velocityField, "DOT_" + velocityField, COL, ts_u_t,
1736 INSERT_VALUES, SCATTER_REVERSE);
1737 CHKERR mField.getInterface<VecManager>()->setOtherGlobalGhostVector(
1738 problemPtr, spatialPositionField, "DOT_" + spatialPositionField, COL,
1739 ts_u_t, INSERT_VALUES, SCATTER_REVERSE);
1740 break;
1741 }
1742 case CTX_TSSETIJACOBIAN: {
1743 snes_B = ts_B;
1744 break;
1745 }
1746 default:
1747 break;
1748 }
1749
1751}
1752
1755 //
1756 // SNES snes;
1757 // CHKERR TSGetSNES(tS,&snes);
1758 // CHKERR SNESSetLagJacobian(snes,jacobianLag);
1760}
1761
1764
1765 Range added_tets;
1767 mField, BLOCKSET | BODYFORCESSET, it)) {
1768 int id = it->getMeshsetId();
1769 EntityHandle meshset = it->getMeshset();
1770 CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET,
1771 setOfBlocks[id].tEts, true);
1772 added_tets.merge(setOfBlocks[id].tEts);
1773 Block_BodyForces mydata;
1774 CHKERR it->getAttributeDataStructure(mydata);
1775 setOfBlocks[id].rho0 = mydata.data.density;
1776 setOfBlocks[id].a0.resize(3);
1777 setOfBlocks[id].a0[0] = mydata.data.acceleration_x;
1778 setOfBlocks[id].a0[1] = mydata.data.acceleration_y;
1779 setOfBlocks[id].a0[2] = mydata.data.acceleration_z;
1780 // std::cerr << setOfBlocks[id].tEts << std::endl;
1781 }
1782
1784 mField, BLOCKSET | MAT_ELASTICSET, it)) {
1785 Mat_Elastic mydata;
1786 CHKERR it->getAttributeDataStructure(mydata);
1787 if (mydata.data.User1 == 0)
1788 continue;
1789 Range tets;
1790 EntityHandle meshset = it->getMeshset();
1791 CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets, true);
1792 tets = subtract(tets, added_tets);
1793 if (tets.empty())
1794 continue;
1795 int id = it->getMeshsetId();
1796 setOfBlocks[-id].tEts = tets;
1797 setOfBlocks[-id].rho0 = mydata.data.User1;
1798 setOfBlocks[-id].a0.resize(3);
1799 setOfBlocks[-id].a0[0] = mydata.data.User2;
1800 setOfBlocks[-id].a0[1] = mydata.data.User3;
1801 setOfBlocks[-id].a0[2] = mydata.data.User4;
1802 // std::cerr << setOfBlocks[id].tEts << std::endl;
1803 }
1804
1806}
1807
1809 MoFEM::Interface &m_field,
1810 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr) {
1812
1813 if (!block_sets_ptr)
1814 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1815 "Pointer to block of sets is null");
1816
1818 m_field, BLOCKSET | BODYFORCESSET, it)) {
1819 Block_BodyForces mydata;
1820 CHKERR it->getAttributeDataStructure(mydata);
1821 int id = it->getMeshsetId();
1822 auto &block_data = (*block_sets_ptr)[id];
1823 EntityHandle meshset = it->getMeshset();
1824 CHKERR m_field.get_moab().get_entities_by_dimension(meshset, 3,
1825 block_data.tEts, true);
1826 block_data.rho0 = mydata.data.density;
1827 block_data.a0.resize(3);
1828 block_data.a0[0] = mydata.data.acceleration_x;
1829 block_data.a0[1] = mydata.data.acceleration_y;
1830 block_data.a0[2] = mydata.data.acceleration_z;
1831 }
1832
1834}
1835
1837 string element_name, string velocity_field_name,
1838 string spatial_position_field_name, string material_position_field_name,
1839 bool ale, BitRefLevel bit) {
1841
1842 //
1843
1844 CHKERR mField.add_finite_element(element_name, MF_ZERO);
1846 velocity_field_name);
1848 velocity_field_name);
1850 velocity_field_name);
1852 element_name, spatial_position_field_name);
1854 element_name, spatial_position_field_name);
1856 element_name, spatial_position_field_name);
1857 if (mField.check_field(material_position_field_name)) {
1858 if (ale) {
1860 element_name, material_position_field_name);
1862 element_name, material_position_field_name);
1864 element_name, "DOT_" + material_position_field_name);
1865 }
1867 element_name, material_position_field_name);
1868 }
1870 element_name, "DOT_" + velocity_field_name);
1872 element_name, "DOT_" + spatial_position_field_name);
1873
1874 Range tets;
1875 if (bit.any()) {
1876 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1877 bit, BitRefLevel().set(), MBTET, tets);
1878 }
1879
1880 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
1881 for (; sit != setOfBlocks.end(); sit++) {
1882 Range add_tets = sit->second.tEts;
1883 if (!tets.empty()) {
1884 add_tets = intersect(add_tets, tets);
1885 }
1887 element_name);
1888 }
1889
1891}
1892
1894 string element_name, string velocity_field_name,
1895 string spatial_position_field_name, string material_position_field_name,
1896 bool ale, BitRefLevel bit) {
1898
1899 //
1900
1901 CHKERR mField.add_finite_element(element_name, MF_ZERO);
1903 velocity_field_name);
1905 velocity_field_name);
1907 velocity_field_name);
1909 element_name, spatial_position_field_name);
1911 element_name, spatial_position_field_name);
1912 if (mField.check_field(material_position_field_name)) {
1913 if (ale) {
1915 element_name, material_position_field_name);
1917 element_name, "DOT_" + material_position_field_name);
1918 }
1920 element_name, material_position_field_name);
1921 }
1923 element_name, "DOT_" + velocity_field_name);
1925 element_name, "DOT_" + spatial_position_field_name);
1926
1927 Range tets;
1928 if (bit.any()) {
1929 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1930 bit, BitRefLevel().set(), MBTET, tets);
1931 }
1932
1933 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
1934 for (; sit != setOfBlocks.end(); sit++) {
1935 Range add_tets = sit->second.tEts;
1936 if (!tets.empty()) {
1937 add_tets = intersect(add_tets, tets);
1938 }
1940 element_name);
1941 }
1942
1944}
1945
1947 string element_name, string velocity_field_name,
1948 string spatial_position_field_name, string material_position_field_name,
1949 bool ale, BitRefLevel bit, Range *intersected) {
1951
1952 //
1953
1954 CHKERR mField.add_finite_element(element_name, MF_ZERO);
1956 velocity_field_name);
1958 velocity_field_name);
1960 element_name, spatial_position_field_name);
1962 element_name, spatial_position_field_name);
1963 if (mField.check_field(material_position_field_name)) {
1964 if (ale) {
1966 element_name, material_position_field_name);
1968 element_name, material_position_field_name);
1970 element_name, "DOT_" + material_position_field_name);
1971 }
1973 element_name, material_position_field_name);
1974 }
1976 element_name, "DOT_" + velocity_field_name);
1978 element_name, "DOT_" + spatial_position_field_name);
1979
1980 Range tets;
1981 if (bit.any()) {
1982 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1983 bit, BitRefLevel().set(), MBTET, tets);
1984 }
1985 if (intersected != NULL) {
1986 if (tets.empty()) {
1987 tets = *intersected;
1988 } else {
1989 tets = intersect(*intersected, tets);
1990 }
1991 }
1992
1993 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
1994 for (; sit != setOfBlocks.end(); sit++) {
1995 Range add_tets = sit->second.tEts;
1996 if (!tets.empty()) {
1997 add_tets = intersect(add_tets, tets);
1998 }
2000 element_name);
2001 }
2002
2004}
2005
2007 string velocity_field_name, string spatial_position_field_name,
2008 string material_position_field_name, bool ale, bool linear) {
2010
2011 commonData.spatialPositions = spatial_position_field_name;
2012 commonData.meshPositions = material_position_field_name;
2013 commonData.spatialVelocities = velocity_field_name;
2014 commonData.lInear = linear;
2015
2016 // Rhs
2017 feMassRhs.getOpPtrVector().push_back(
2018 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2019 feMassRhs.getOpPtrVector().push_back(
2020 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2021 feMassRhs.getOpPtrVector().push_back(
2022 new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2024 "DOT_" + spatial_position_field_name, commonData));
2025 if (mField.check_field(material_position_field_name)) {
2027 material_position_field_name, commonData));
2028 if (ale) {
2030 "DOT_" + material_position_field_name, commonData));
2031 } else {
2032 feMassRhs.meshPositionsFieldName = material_position_field_name;
2033 }
2034 }
2035 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
2036 for (; sit != setOfBlocks.end(); sit++) {
2037 feMassRhs.getOpPtrVector().push_back(
2038 new OpMassJacobian(spatial_position_field_name, sit->second, commonData,
2039 methodsOp, tAg, false));
2040 feMassRhs.getOpPtrVector().push_back(
2041 new OpMassRhs(spatial_position_field_name, sit->second, commonData));
2042 }
2043
2044 // Lhs
2045 feMassLhs.getOpPtrVector().push_back(
2046 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2047 feMassLhs.getOpPtrVector().push_back(
2048 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2049 feMassLhs.getOpPtrVector().push_back(
2050 new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2052 "DOT_" + spatial_position_field_name, commonData));
2053 if (mField.check_field(material_position_field_name)) {
2055 material_position_field_name, commonData));
2056 if (ale) {
2058 "DOT_" + material_position_field_name, commonData));
2059 } else {
2060 feMassLhs.meshPositionsFieldName = material_position_field_name;
2061 }
2062 }
2063 sit = setOfBlocks.begin();
2064 for (; sit != setOfBlocks.end(); sit++) {
2065 feMassLhs.getOpPtrVector().push_back(
2066 new OpMassJacobian(spatial_position_field_name, sit->second, commonData,
2067 methodsOp, tAg, true));
2068 feMassLhs.getOpPtrVector().push_back(
2069 new OpMassLhs_dM_dv(spatial_position_field_name, velocity_field_name,
2070 sit->second, commonData));
2072 spatial_position_field_name, spatial_position_field_name, sit->second,
2073 commonData));
2074 if (mField.check_field(material_position_field_name)) {
2075 if (ale) {
2077 spatial_position_field_name, material_position_field_name,
2078 sit->second, commonData));
2079 } else {
2080 feMassLhs.meshPositionsFieldName = material_position_field_name;
2081 }
2082 }
2083 }
2084
2085 // Energy
2086 feEnergy.getOpPtrVector().push_back(
2087 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2088 feEnergy.getOpPtrVector().push_back(
2089 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2090 if (mField.check_field(material_position_field_name)) {
2092 material_position_field_name, commonData));
2093 feEnergy.meshPositionsFieldName = material_position_field_name;
2094 }
2095 sit = setOfBlocks.begin();
2096 for (; sit != setOfBlocks.end(); sit++) {
2097 feEnergy.getOpPtrVector().push_back(new OpEnergy(
2098 spatial_position_field_name, sit->second, commonData, feEnergy.V));
2099 }
2100
2102}
2103
2105 string velocity_field_name, string spatial_position_field_name,
2106 string material_position_field_name, bool ale) {
2108
2109 commonData.spatialPositions = spatial_position_field_name;
2110 commonData.meshPositions = material_position_field_name;
2111 commonData.spatialVelocities = velocity_field_name;
2112
2113 // Rhs
2114 feVelRhs.getOpPtrVector().push_back(
2115 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2116 feVelRhs.getOpPtrVector().push_back(
2117 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2118 feVelRhs.getOpPtrVector().push_back(
2119 new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2120 if (mField.check_field(material_position_field_name)) {
2122 "DOT_" + spatial_position_field_name, commonData));
2124 material_position_field_name, commonData));
2125 if (ale) {
2127 "DOT_" + material_position_field_name, commonData));
2128 } else {
2129 feVelRhs.meshPositionsFieldName = material_position_field_name;
2130 }
2131 }
2132 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
2133 for (; sit != setOfBlocks.end(); sit++) {
2135 velocity_field_name, sit->second, commonData, tAg, false));
2136 feVelRhs.getOpPtrVector().push_back(
2137 new OpVelocityRhs(velocity_field_name, sit->second, commonData));
2138 }
2139
2140 // Lhs
2141 feVelLhs.getOpPtrVector().push_back(
2142 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2143 feVelLhs.getOpPtrVector().push_back(
2144 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2145 feVelLhs.getOpPtrVector().push_back(
2146 new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2147 if (mField.check_field(material_position_field_name)) {
2149 "DOT_" + spatial_position_field_name, commonData));
2151 material_position_field_name, commonData));
2152 if (ale) {
2154 "DOT_" + material_position_field_name, commonData));
2155 } else {
2156 feVelLhs.meshPositionsFieldName = material_position_field_name;
2157 }
2158 }
2159 sit = setOfBlocks.begin();
2160 for (; sit != setOfBlocks.end(); sit++) {
2162 velocity_field_name, sit->second, commonData, tAg));
2164 velocity_field_name, velocity_field_name, sit->second, commonData));
2166 velocity_field_name, spatial_position_field_name, sit->second,
2167 commonData));
2168 if (mField.check_field(material_position_field_name)) {
2169 if (ale) {
2171 velocity_field_name, material_position_field_name, sit->second,
2172 commonData));
2173 } else {
2174 feVelLhs.meshPositionsFieldName = material_position_field_name;
2175 }
2176 }
2177 }
2178
2180}
2181
2183 string velocity_field_name, string spatial_position_field_name,
2184 string material_position_field_name, Range *forces_on_entities_ptr) {
2186
2187 commonData.spatialPositions = spatial_position_field_name;
2188 commonData.meshPositions = material_position_field_name;
2189 commonData.spatialVelocities = velocity_field_name;
2190
2191 // Rhs
2192 feTRhs.getOpPtrVector().push_back(
2193 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2194 feTRhs.getOpPtrVector().push_back(
2195 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2196 feTRhs.getOpPtrVector().push_back(
2197 new OpGetCommonDataAtGaussPts(material_position_field_name, commonData));
2198 feTRhs.getOpPtrVector().push_back(
2199 new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2200
2201 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
2202 for (; sit != setOfBlocks.end(); sit++) {
2203 feTRhs.getOpPtrVector().push_back(
2205 material_position_field_name, sit->second, commonData, tAg, false));
2207 material_position_field_name, sit->second, commonData,
2208 forces_on_entities_ptr));
2209 }
2210
2211 // Lhs
2212 feTLhs.getOpPtrVector().push_back(
2213 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2214 feTLhs.getOpPtrVector().push_back(
2215 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2216 feTLhs.getOpPtrVector().push_back(
2217 new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2218 if (mField.check_field(material_position_field_name)) {
2220 material_position_field_name, commonData));
2221 }
2222 sit = setOfBlocks.begin();
2223 for (; sit != setOfBlocks.end(); sit++) {
2224 feTLhs.getOpPtrVector().push_back(
2226 material_position_field_name, sit->second, commonData, tAg));
2227 feTLhs.getOpPtrVector().push_back(
2229 material_position_field_name, velocity_field_name, sit->second,
2230 commonData, forces_on_entities_ptr));
2231 feTLhs.getOpPtrVector().push_back(
2233 material_position_field_name, spatial_position_field_name,
2234 sit->second, commonData, forces_on_entities_ptr));
2235 feTLhs.getOpPtrVector().push_back(
2237 material_position_field_name, material_position_field_name,
2238 sit->second, commonData, forces_on_entities_ptr));
2239 }
2240
2242}
2243
2245 string velocity_field_name, string spatial_position_field_name,
2246 string material_position_field_name, bool linear) {
2248
2249 commonData.spatialPositions = spatial_position_field_name;
2250 commonData.meshPositions = material_position_field_name;
2251 commonData.spatialVelocities = velocity_field_name;
2252 commonData.lInear = linear;
2253
2254 // Rhs
2255 feMassRhs.getOpPtrVector().push_back(
2256 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2257 feMassRhs.getOpPtrVector().push_back(
2258 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2259 feMassRhs.getOpPtrVector().push_back(
2260 new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2261 if (mField.check_field(material_position_field_name)) {
2263 material_position_field_name, commonData));
2264 feMassRhs.meshPositionsFieldName = material_position_field_name;
2265 }
2266 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
2267 for (; sit != setOfBlocks.end(); sit++) {
2268 feMassRhs.getOpPtrVector().push_back(
2269 new OpMassJacobian(spatial_position_field_name, sit->second, commonData,
2270 methodsOp, tAg, false));
2271 feMassRhs.getOpPtrVector().push_back(
2272 new OpMassRhs(spatial_position_field_name, sit->second, commonData));
2273 }
2274
2275 // Lhs
2276 feMassLhs.getOpPtrVector().push_back(
2277 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2278 feMassLhs.getOpPtrVector().push_back(
2279 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2280 feMassLhs.getOpPtrVector().push_back(
2281 new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2282 if (mField.check_field(material_position_field_name)) {
2284 material_position_field_name, commonData));
2285 feMassLhs.meshPositionsFieldName = material_position_field_name;
2286 }
2287 sit = setOfBlocks.begin();
2288 for (; sit != setOfBlocks.end(); sit++) {
2289 feMassLhs.getOpPtrVector().push_back(
2290 new OpMassJacobian(spatial_position_field_name, sit->second, commonData,
2291 methodsOp, tAg, true));
2293 spatial_position_field_name, spatial_position_field_name, sit->second,
2294 commonData));
2295 if (mField.check_field(material_position_field_name)) {
2296 feMassLhs.meshPositionsFieldName = material_position_field_name;
2297 }
2298 }
2299
2300 // Aux Lhs
2301 feMassAuxLhs.getOpPtrVector().push_back(
2302 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2303 feMassAuxLhs.getOpPtrVector().push_back(
2304 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2305 feMassAuxLhs.getOpPtrVector().push_back(
2306 new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2307 if (mField.check_field(material_position_field_name)) {
2309 material_position_field_name, commonData));
2310 feMassAuxLhs.meshPositionsFieldName = material_position_field_name;
2311 }
2312 sit = setOfBlocks.begin();
2313 for (; sit != setOfBlocks.end(); sit++) {
2314 feMassAuxLhs.getOpPtrVector().push_back(
2315 new OpMassJacobian(spatial_position_field_name, sit->second, commonData,
2316 methodsOp, tAg, true));
2318 spatial_position_field_name, spatial_position_field_name, sit->second,
2319 commonData));
2320 if (mField.check_field(material_position_field_name)) {
2321 feMassAuxLhs.meshPositionsFieldName = material_position_field_name;
2322 }
2323 }
2324
2325 // Energy E=0.5*rho*v*v
2326 feEnergy.getOpPtrVector().push_back(
2327 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2328 feEnergy.getOpPtrVector().push_back(
2329 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2330 if (mField.check_field(material_position_field_name)) {
2332 material_position_field_name, commonData));
2333 feEnergy.meshPositionsFieldName = material_position_field_name;
2334 }
2335 sit = setOfBlocks.begin();
2336 for (; sit != setOfBlocks.end(); sit++) {
2337 feEnergy.getOpPtrVector().push_back(new OpEnergy(
2338 spatial_position_field_name, sit->second, commonData, feEnergy.V));
2339 }
2340
2342}
2343
2346 if (iNitialized) {
2347
2348 CHKERR dEstroy();
2349 CHKERRABORT(PETSC_COMM_WORLD, ierr);
2350 }
2351}
2352
2355 if (!iNitialized) {
2356
2357#if PETSC_VERSION_GE(3, 5, 3)
2358 CHKERR MatCreateVecs(K, &u, &Ku);
2359 CHKERR MatCreateVecs(M, &v, &Mv);
2360#else
2361 CHKERR MatGetVecs(K, &u, &Ku);
2362 CHKERR MatGetVecs(M, &v, &Mv);
2363#endif
2364 CHKERR MatDuplicate(K, MAT_SHARE_NONZERO_PATTERN, &barK);
2365 iNitialized = true;
2366 }
2368}
2369
2372 if (iNitialized) {
2373
2374 CHKERR VecDestroy(&u);
2375 CHKERR VecDestroy(&Ku);
2376 CHKERR VecDestroy(&v);
2377 CHKERR VecDestroy(&Mv);
2378 CHKERR MatDestroy(&barK);
2379 iNitialized = false;
2380 }
2382}
2383
2386
2387 if (!initPC) {
2388 MPI_Comm comm;
2389 CHKERR PetscObjectGetComm((PetscObject)shellMat, &comm);
2390 CHKERR PCCreate(comm, &pC);
2391 initPC = true;
2392 }
2394}
2395
2398
2399 if (initPC) {
2400 CHKERR PCDestroy(&pC);
2401 initPC = false;
2402 }
2404}
2405
2409
2412
2413 if (ts_ctx != CTX_TSSETIFUNCTION) {
2414 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
2415 "It is used to residual of velocities");
2416 }
2417 if (!shellMatCtx->iNitialized) {
2418 CHKERR shellMatCtx->iNit();
2419 }
2420 // Note velocities calculate from displacements are stroed in shellMatCtx->u
2421 CHKERR VecScatterBegin(shellMatCtx->scatterU, ts_u_t, shellMatCtx->u,
2422 INSERT_VALUES, SCATTER_FORWARD);
2423 CHKERR VecScatterEnd(shellMatCtx->scatterU, ts_u_t, shellMatCtx->u,
2424 INSERT_VALUES, SCATTER_FORWARD);
2425 CHKERR VecScatterBegin(shellMatCtx->scatterV, ts_u, shellMatCtx->v,
2426 INSERT_VALUES, SCATTER_FORWARD);
2427 CHKERR VecScatterEnd(shellMatCtx->scatterV, ts_u, shellMatCtx->v,
2428 INSERT_VALUES, SCATTER_FORWARD);
2429 CHKERR VecAXPY(shellMatCtx->v, -1, shellMatCtx->u);
2430 CHKERR VecScatterBegin(shellMatCtx->scatterV, shellMatCtx->v, ts_F,
2431 ADD_VALUES, SCATTER_REVERSE);
2432 CHKERR VecScatterEnd(shellMatCtx->scatterV, shellMatCtx->v, ts_F, ADD_VALUES,
2433 SCATTER_REVERSE);
2434 // VecView(shellMatCtx->v,PETSC_VIEWER_STDOUT_WORLD);
2435
2437}
2438
2439#ifdef __DIRICHLET_HPP__
2440
2441ConvectiveMassElement::ShellMatrixElement::ShellMatrixElement(
2442 MoFEM::Interface &m_field)
2443 : mField(m_field) {}
2444
2445MoFEMErrorCode ConvectiveMassElement::ShellMatrixElement::preProcess() {
2447
2448 if (ts_ctx != CTX_TSSETIJACOBIAN) {
2449 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
2450 "It is used to calculate shell matrix only");
2451 }
2452
2453 shellMatCtx->ts_a = ts_a;
2454 DirichletBcPtr->copyTs(*((TSMethod *)this)); // copy context for TSMethod
2455
2456 DirichletBcPtr->dIag = 1;
2457 DirichletBcPtr->ts_B = shellMatCtx->K;
2458 CHKERR MatZeroEntries(shellMatCtx->K);
2459 CHKERR mField.problem_basic_method_preProcess(problemName, *DirichletBcPtr);
2460 LoopsToDoType::iterator itk = loopK.begin();
2461 for (; itk != loopK.end(); itk++) {
2462 itk->second->copyTs(*((TSMethod *)this));
2463 itk->second->ts_B = shellMatCtx->K;
2464 CHKERR mField.loop_finite_elements(problemName, itk->first, *itk->second);
2465 }
2466 LoopsToDoType::iterator itam = loopAuxM.begin();
2467 for (; itam != loopAuxM.end(); itam++) {
2468 itam->second->copyTs(*((TSMethod *)this));
2469 itam->second->ts_B = shellMatCtx->K;
2470 CHKERR mField.loop_finite_elements(problemName, itam->first, *itam->second);
2471 }
2472 CHKERR mField.problem_basic_method_postProcess(problemName, *DirichletBcPtr);
2473 CHKERR MatAssemblyBegin(shellMatCtx->K, MAT_FINAL_ASSEMBLY);
2474 CHKERR MatAssemblyEnd(shellMatCtx->K, MAT_FINAL_ASSEMBLY);
2475
2476 DirichletBcPtr->dIag = 0;
2477 DirichletBcPtr->ts_B = shellMatCtx->M;
2478 CHKERR MatZeroEntries(shellMatCtx->M);
2479 // CHKERR mField.problem_basic_method_preProcess(problemName,*DirichletBcPtr);
2480 LoopsToDoType::iterator itm = loopM.begin();
2481 for (; itm != loopM.end(); itm++) {
2482 itm->second->copyTs(*((TSMethod *)this));
2483 itm->second->ts_B = shellMatCtx->M;
2484 CHKERR mField.loop_finite_elements(problemName, itm->first, *itm->second);
2485 }
2486 CHKERR mField.problem_basic_method_postProcess(problemName, *DirichletBcPtr);
2487 CHKERR MatAssemblyBegin(shellMatCtx->M, MAT_FINAL_ASSEMBLY);
2488 CHKERR MatAssemblyEnd(shellMatCtx->M, MAT_FINAL_ASSEMBLY);
2489
2490 // barK
2491 CHKERR MatZeroEntries(shellMatCtx->barK);
2492 CHKERR MatCopy(shellMatCtx->K, shellMatCtx->barK, SAME_NONZERO_PATTERN);
2493 CHKERR MatAXPY(shellMatCtx->barK, ts_a, shellMatCtx->M, SAME_NONZERO_PATTERN);
2494 CHKERR MatAssemblyBegin(shellMatCtx->barK, MAT_FINAL_ASSEMBLY);
2495 CHKERR MatAssemblyEnd(shellMatCtx->barK, MAT_FINAL_ASSEMBLY);
2496
2497 // Matrix View
2498 // MatView(shellMatCtx->barK,PETSC_VIEWER_DRAW_WORLD);//PETSC_VIEWER_STDOUT_WORLD);
2499 // std::string wait;
2500 // std::cin >> wait;
2501
2503}
2504
2505#endif //__DIRICHLET_HPP__
Operators and data structures for mass and convective mass element.
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr double a
@ COL
@ MF_ZERO
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ BODYFORCESSET
block name is "BODY_FORCES"
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
@ BLOCKSET
@ 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()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
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_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
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
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
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.
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.
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
auto getFTensor2FromArray3by3(ublas::matrix< T, L, A > &data, const FTensor::Number< S > &, const size_t rr, const size_t cc=0)
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
double h
double rho
Definition plastic.cpp:144
double H
Hardening.
Definition plastic.cpp:128
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
MoFEMErrorCode postProcess()
function is run at the end of loop
int getRule(int order)
it is used to calculate nb. of Gauss integration points
MoFEMErrorCode preProcess()
function is run at the beginning of loop
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 postProcess()
function is run at the end of loop
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 field data coefficients.
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
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Elastic material data structure.
intrusive_ptr for managing petsc objects
data structure for TS (time stepping) context
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Vector manager is used to create vectors \mofem_vectors.