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);
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)
constexpr IntegrationType G
Definition level_set.cpp:33
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.
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.
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
constexpr AssemblyType A
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()
Post-processing function executed at loop completion.
int getRule(int order)
it is used to calculate nb. of Gauss integration points
MoFEMErrorCode preProcess()
Pre-processing function executed at loop initialization.
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()
Post-processing function executed at loop completion.
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 DOF values on entity.
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 DOF data structures (const version)
const VectorInt & getIndices() const
Get global indices of degrees of freedom on entity.
structure to get information from 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.