v0.13.1
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
8
9#include <MoFEM.hpp>
10using namespace MoFEM;
11
13
14#include <adolc/adolc.h>
16#include <DirichletBC.hpp>
19
20#ifndef WITH_ADOL_C
21#error "MoFEM need to be compiled with ADOL-C"
22#endif
23
25 : VolumeElementForcesAndSourcesCore(m_field), A(PETSC_NULL), F(PETSC_NULL) {
26 meshPositionsFieldName = "NoNE";
27
28 auto create_vec = [&]() {
29 constexpr int ghosts[] = {0};
30 if (mField.get_comm_rank() == 0) {
32 } else {
33 return createSmartVectorMPI(mField.get_comm(), 0, 1);
34 }
35 };
36
37 V = create_vec();
38}
39
41
44
45 CHKERR VolumeElementForcesAndSourcesCore::preProcess();
46
47 switch (ts_ctx) {
48 case CTX_TSNONE:
49 CHKERR VecZeroEntries(V);
50 break;
51 default:
52 break;
53 }
54
56}
57
60
61 CHKERR VolumeElementForcesAndSourcesCore::postProcess();
62
63 const double *array;
64 switch (ts_ctx) {
65 case CTX_TSNONE:
66 CHKERR VecAssemblyBegin(V);
67 CHKERR VecAssemblyEnd(V);
68 CHKERR VecSum(V, &eNergy);
69 break;
70 default:
71 break;
72 }
73
75}
76
78 short int tag)
79 : feMassRhs(m_field), feMassLhs(m_field), feMassAuxLhs(m_field),
80 feVelRhs(m_field), feVelLhs(m_field), feTRhs(m_field), feTLhs(m_field),
81 feEnergy(m_field), mField(m_field), tAg(tag) {}
82
84 const std::string field_name,
85 std::vector<VectorDouble> &values_at_gauss_pts,
86 std::vector<MatrixDouble> &gardient_at_gauss_pts)
89 valuesAtGaussPts(values_at_gauss_pts),
90 gradientAtGaussPts(gardient_at_gauss_pts), zeroAtType(MBVERTEX) {}
91
95
96 int nb_dofs = data.getFieldData().size();
97 if (nb_dofs == 0) {
99 }
100 int nb_gauss_pts = data.getN().size1();
101 int nb_base_functions = data.getN().size2();
102
103 // initialize
104 // VectorDouble& values = data.getFieldData();
105 valuesAtGaussPts.resize(nb_gauss_pts);
106 gradientAtGaussPts.resize(nb_gauss_pts);
107 for (int gg = 0; gg < nb_gauss_pts; gg++) {
108 valuesAtGaussPts[gg].resize(3);
109 gradientAtGaussPts[gg].resize(3, 3);
110 }
111
112 if (type == zeroAtType) {
113 for (int gg = 0; gg < nb_gauss_pts; gg++) {
114 valuesAtGaussPts[gg].clear();
115 gradientAtGaussPts[gg].clear();
116 }
117 }
118
119 auto base_function = data.getFTensor0N();
120 auto diff_base_functions = data.getFTensor1DiffN<3>();
123
124 for (int gg = 0; gg != nb_gauss_pts; gg++) {
125 auto field_data = data.getFTensor1FieldData<3>();
126 FTensor::Tensor1<double *, 3> values(&valuesAtGaussPts[gg][0],
127 &valuesAtGaussPts[gg][1],
128 &valuesAtGaussPts[gg][2]);
130 &gradientAtGaussPts[gg](0, 0), &gradientAtGaussPts[gg](0, 1),
131 &gradientAtGaussPts[gg](0, 2), &gradientAtGaussPts[gg](1, 0),
132 &gradientAtGaussPts[gg](1, 1), &gradientAtGaussPts[gg](1, 2),
133 &gradientAtGaussPts[gg](2, 0), &gradientAtGaussPts[gg](2, 1),
134 &gradientAtGaussPts[gg](2, 2));
135 int bb = 0;
136 for (; bb != nb_dofs / 3; bb++) {
137 values(i) += base_function * field_data(i);
138 gradient(i, j) += field_data(i) * diff_base_functions(j);
139 ++diff_base_functions;
140 ++base_function;
141 ++field_data;
142 }
143 for (; bb != nb_base_functions; bb++) {
144 ++diff_base_functions;
145 ++base_function;
146 }
147 }
149}
150
152 const std::string field_name, CommonData &common_data)
153 : OpGetDataAtGaussPts(field_name, common_data.dataAtGaussPts[field_name],
154 common_data.gradAtGaussPts[field_name]) {}
155
157 const std::string field_name, BlockData &data, CommonData &common_data,
158 boost::ptr_vector<MethodForForceScaling> &methods_op, int tag,
159 bool jacobian)
162 dAta(data), commonData(common_data), tAg(tag), jAcobian(jacobian),
163 lInear(commonData.lInear), fieldDisp(false), methodsOp(methods_op) {}
164
166 int row_side, EntityType row_type,
167 EntitiesFieldData::EntData &row_data) {
169
170 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
171 dAta.tEts.end()) {
173 }
174
175 // do it only once, no need to repeat this for edges,faces or tets
176 if (row_type != MBVERTEX)
178
179 int nb_dofs = row_data.getIndices().size();
180 if (nb_dofs == 0)
182
183 {
184
185 if (a.size() != 3) {
186 a.resize(3, false);
187 dot_W.resize(3, false);
188 a_res.resize(3, false);
189 g.resize(3, 3, false);
190 G.resize(3, 3, false);
191 h.resize(3, 3, false);
192 H.resize(3, 3, false);
193 invH.resize(3, 3, false);
194 F.resize(3, 3, false);
195 }
196
197 dot_W.clear();
198 H.clear();
199 invH.clear();
200 for (int dd = 0; dd < 3; dd++) {
201 H(dd, dd) = 1;
202 invH(dd, dd) = 1;
203 }
204
205 int nb_gauss_pts = row_data.getN().size1();
206 commonData.valMass.resize(nb_gauss_pts);
207 commonData.jacMassRowPtr.resize(nb_gauss_pts);
208 commonData.jacMass.resize(nb_gauss_pts);
209
210 const std::vector<VectorDouble> &dot_spacial_vel =
212
213 const std::vector<MatrixDouble> &spatial_positions_grad =
215
216 const std::vector<MatrixDouble> &spatial_velocities_grad =
218
219 const std::vector<VectorDouble> &meshpos_vel =
221
222 const std::vector<MatrixDouble> &mesh_positions_gradient =
224
225 int nb_active_vars = 0;
226 for (int gg = 0; gg < nb_gauss_pts; gg++) {
227
228 if (gg == 0) {
229
230 trace_on(tAg);
231
232 for (int nn1 = 0; nn1 < 3; nn1++) { // 0
233 // commonData.dataAtGaussPts["DOT_"+commonData.spatialVelocities]
234 a[nn1] <<= dot_spacial_vel[gg][nn1];
235 nb_active_vars++;
236 }
237 for (int nn1 = 0; nn1 < 3; nn1++) { // 3
238 for (int nn2 = 0; nn2 < 3; nn2++) {
239 // commonData.gradAtGaussPts[commonData.spatialPositions][gg]
240 h(nn1, nn2) <<= spatial_positions_grad[gg](nn1, nn2);
241 if (fieldDisp) {
242 if (nn1 == nn2) {
243 h(nn1, nn2) += 1;
244 }
245 }
246 nb_active_vars++;
247 }
248 }
250 .size() > 0) {
251 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+9=12
252 for (int nn2 = 0; nn2 < 3; nn2++) {
253 // commonData.gradAtGaussPts[commonData.spatialVelocities]
254 g(nn1, nn2) <<= spatial_velocities_grad[gg](nn1, nn2);
255 nb_active_vars++;
256 }
257 }
258 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+9+9=21
259 // commonData.dataAtGaussPts["DOT_"+commonData.meshPositions]
260 dot_W(nn1) <<= meshpos_vel[gg][nn1];
261 nb_active_vars++;
262 }
263 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+9+9+3=24
264 for (int nn2 = 0; nn2 < 3; nn2++) {
265 // commonData.gradAtGaussPts[commonData.meshPositions][gg]
266 H(nn1, nn2) <<= mesh_positions_gradient[gg](nn1, nn2);
267 nb_active_vars++;
268 }
269 }
270 }
271
272 auto a0 = dAta.a0;
274
275 auto t_a_res =
276 FTensor::Tensor1<adouble *, 3>{&a_res[0], &a_res[1], &a_res[2]};
277 auto t_a = FTensor::Tensor1<adouble *, 3>{&a[0], &a[1], &a[2]};
278 auto t_a0 = FTensor::Tensor1<double *, 3>{&a0[0], &a0[1], &a0[2]};
279 auto t_dotW =
280 FTensor::Tensor1<adouble *, 3>{&dot_W[0], &dot_W[1], &dot_W[2]};
283 auto t_invH = getFTensor2FromArray3by3(invH, FTensor::Number<0>(), 0);
286
287 const double rho0 = dAta.rho0;
288
290 CHKERR invertTensor3by3(H, detH, invH);
291
292 t_G(i, j) = t_g(i, k) * t_invH(k, j);
293 t_a_res(i) = t_a(i) - t_a0(i) + t_G(i, j) * t_dotW(j);
294
295 if (!lInear) {
296
297 t_F(i,j) = t_h(i,k)*t_invH(k,j);
298 t_a_res(i) *= rho0 * detH;
299 t_a_res(i) *= determinantTensor3by3(t_F);
300
301 } else {
302
303 t_a_res(i) *= rho0 * detH;
304
305 }
306
307 // dependant
309 res.resize(3);
310 for (int rr = 0; rr < 3; rr++) {
311 a_res[rr] >>= res[rr];
312 }
313
314 trace_off();
315 }
316
317 active.resize(nb_active_vars);
318 int aa = 0;
319 for (int nn1 = 0; nn1 < 3; nn1++) { // 0
320 active[aa++] = dot_spacial_vel[gg][nn1];
321 }
322 for (int nn1 = 0; nn1 < 3; nn1++) { // 3
323 for (int nn2 = 0; nn2 < 3; nn2++) {
324 if (fieldDisp && nn1 == nn2) {
325 active[aa++] = spatial_positions_grad[gg](nn1, nn2) + 1;
326 } else {
327 active[aa++] = spatial_positions_grad[gg](nn1, nn2);
328 }
329 }
330 }
332 0) {
333 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+9=12
334 for (int nn2 = 0; nn2 < 3; nn2++) {
335 active[aa++] = spatial_velocities_grad[gg](nn1, nn2);
336 }
337 }
338 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+9+9=21
339 active[aa++] = meshpos_vel[gg][nn1];
340 }
341 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+9+9+3=24
342 for (int nn2 = 0; nn2 < 3; nn2++) {
343 active[aa++] = mesh_positions_gradient[gg](nn1, nn2);
344 }
345 }
346 }
347
348 if (!jAcobian) {
350 if (gg > 0) {
351 res.resize(3);
352 int r;
353 r = ::function(tAg, 3, nb_active_vars, &active[0], &res[0]);
354 if (r != 3) { // function is locally analytic
355 SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
356 "ADOL-C function evaluation with error r = %d", r);
357 }
358 }
359 double val = getVolume() * getGaussPts()(3, gg);
360 res *= val;
361 } else {
362 commonData.jacMassRowPtr[gg].resize(3);
363 commonData.jacMass[gg].resize(3, nb_active_vars);
364 for (int nn1 = 0; nn1 < 3; nn1++) {
365 (commonData.jacMassRowPtr[gg])[nn1] =
366 &(commonData.jacMass[gg](nn1, 0));
367 }
368 int r;
369 r = jacobian(tAg, 3, nb_active_vars, &active[0],
370 &(commonData.jacMassRowPtr[gg])[0]);
371 if (r != 3) {
372 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
373 "ADOL-C function evaluation with error");
374 }
375 double val = getVolume() * getGaussPts()(3, gg);
376 commonData.jacMass[gg] *= val;
377 }
378 }
379 }
380
382}
383
385 BlockData &data,
386 CommonData &common_data)
389 dAta(data), commonData(common_data) {}
390
392 int row_side, EntityType row_type,
393 EntitiesFieldData::EntData &row_data) {
395
396 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
397 dAta.tEts.end()) {
399 }
400 if (row_data.getIndices().size() == 0)
402 int nb_dofs = row_data.getIndices().size();
403
404 auto base = row_data.getFTensor0N();
405 int nb_base_functions = row_data.getN().size2();
406
407 {
408
409 nf.resize(nb_dofs);
410 nf.clear();
411
413
414 for (unsigned int gg = 0; gg < row_data.getN().size1(); gg++) {
415 FTensor::Tensor1<double *, 3> t_nf(&nf[0], &nf[1], &nf[2], 3);
417 &commonData.valMass[gg][1],
418 &commonData.valMass[gg][2]);
419 int dd = 0;
420 for (; dd < nb_dofs / 3; dd++) {
421 t_nf(i) += base * res(i);
422 ++base;
423 ++t_nf;
424 }
425 for (; dd != nb_base_functions; dd++) {
426 ++base;
427 }
428 }
429
430 if ((unsigned int)nb_dofs > 3 * row_data.getN().size2()) {
431 SETERRQ(PETSC_COMM_SELF, 1, "data inconsistency");
432 }
433 CHKERR VecSetValues(getFEMethod()->ts_F, nb_dofs, &row_data.getIndices()[0],
434 &nf[0], ADD_VALUES);
435 }
436
438}
439
441 const std::string vel_field, const std::string field_name, BlockData &data,
442 CommonData &common_data, Range *forcesonlyonentities_ptr)
444 vel_field, field_name,
446 dAta(data), commonData(common_data) {
447 sYmm = false;
448 if (forcesonlyonentities_ptr != NULL) {
449 forcesOnlyOnEntities = *forcesonlyonentities_ptr;
450 }
451}
452
454 EntitiesFieldData::EntData &col_data, int gg) {
456 int nb_col = col_data.getIndices().size();
457 jac.clear();
458 if (!nb_col)
463 FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
464 &jac(1, 0), &jac(1, 1), &jac(1, 2),
465 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
467 &commonData.jacMass[gg](0, 0), &commonData.jacMass[gg](0, 1),
468 &commonData.jacMass[gg](0, 2), &commonData.jacMass[gg](1, 0),
469 &commonData.jacMass[gg](1, 1), &commonData.jacMass[gg](1, 2),
470 &commonData.jacMass[gg](2, 0), &commonData.jacMass[gg](2, 1),
471 &commonData.jacMass[gg](2, 2));
472 double *base_ptr = const_cast<double *>(&col_data.getN(gg)[0]);
473 FTensor::Tensor0<double *> base(base_ptr, 1);
474 if (commonData.dataAtGaussPts["DOT_" + commonData.meshPositions].size() ==
475 0) {
476 for (int dd = 0; dd < nb_col / 3; dd++) {
477 t_jac(i, j) += t_mass1(i, j) * base * getFEMethod()->ts_a;
478 ++base;
479 ++t_jac;
480 }
481 } else {
482 const int s = 3 + 9;
484 // T* d000, T* d001, T* d002,
485 // T* d010, T* d011, T* d012,
486 // T* d020, T* d021, T* d022,
487 // T* d100, T* d101, T* d102,
488 // T* d110, T* d111, T* d112,
489 // T* d120, T* d121, T* d122,
490 // T* d200, T* d201, T* d202,
491 // T* d210, T* d211, T* d212,
492 // T* d220, T* d221, T* d222,
493 &commonData.jacMass[gg](0, s + 0), &commonData.jacMass[gg](0, s + 1),
494 &commonData.jacMass[gg](0, s + 2), &commonData.jacMass[gg](0, s + 3),
495 &commonData.jacMass[gg](0, s + 4), &commonData.jacMass[gg](0, s + 5),
496 &commonData.jacMass[gg](0, s + 6), &commonData.jacMass[gg](0, s + 7),
497 &commonData.jacMass[gg](0, s + 8), &commonData.jacMass[gg](1, s + 0),
498 &commonData.jacMass[gg](1, s + 1), &commonData.jacMass[gg](1, s + 2),
499 &commonData.jacMass[gg](1, s + 3), &commonData.jacMass[gg](1, s + 4),
500 &commonData.jacMass[gg](1, s + 5), &commonData.jacMass[gg](1, s + 6),
501 &commonData.jacMass[gg](1, s + 7), &commonData.jacMass[gg](1, s + 8),
502 &commonData.jacMass[gg](2, s + 0), &commonData.jacMass[gg](2, s + 1),
503 &commonData.jacMass[gg](2, s + 2), &commonData.jacMass[gg](2, s + 3),
504 &commonData.jacMass[gg](2, s + 4), &commonData.jacMass[gg](2, s + 5),
505 &commonData.jacMass[gg](2, s + 6), &commonData.jacMass[gg](2, s + 7),
506 &commonData.jacMass[gg](2, s + 8));
507
508 double *diff_ptr =
509 const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
510 FTensor::Tensor1<double *, 3> diff(diff_ptr, &diff_ptr[1], &diff_ptr[2], 3);
511 for (int dd = 0; dd < nb_col / 3; dd++) {
512 t_jac(i, j) += t_mass1(i, j) * base * getFEMethod()->ts_a;
513 t_jac(i, j) += t_mass3(i, j, k) * diff(k);
514 ++base;
515 ++diff;
516 ++t_jac;
517 }
518 }
520}
521
523 int row_side, int col_side, EntityType row_type, EntityType col_type,
525 EntitiesFieldData::EntData &col_data) {
527
528 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
529 dAta.tEts.end()) {
531 }
532
533 int nb_row = row_data.getIndices().size();
534 int nb_col = col_data.getIndices().size();
535 if (nb_row == 0)
537 if (nb_col == 0)
539
540 auto base = row_data.getFTensor0N();
541 int nb_base_functions = row_data.getN().size2();
542
543 {
544
545 k.resize(nb_row, nb_col);
546 k.clear();
547 jac.resize(3, nb_col);
548
549 for (unsigned int gg = 0; gg < row_data.getN().size1(); gg++) {
550
551 try {
552 CHKERR getJac(col_data, gg);
553 } catch (const std::exception &ex) {
554 std::ostringstream ss;
555 ss << "throw in method: " << ex.what() << std::endl;
556 SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
557 }
558
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) {
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));
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
738 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,
799 EntitiesFieldData::EntData &row_data) {
801
802 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
803 dAta.tEts.end()) {
805 }
806
807 // do it only once, no need to repeat this for edges,faces or tets
808 if (row_type != MBVERTEX)
810
811 int nb_dofs = row_data.getIndices().size();
812 if (nb_dofs == 0)
814
815 {
816
817 v.resize(3);
818 dot_w.resize(3);
819 h.resize(3, 3);
820 h.clear();
821 F.resize(3, 3);
822 dot_W.resize(3);
823 dot_W.clear();
824 H.resize(3, 3);
825 H.clear();
826 invH.resize(3, 3);
827 invH.clear();
828 dot_u.resize(3);
829 for (int dd = 0; dd < 3; dd++) {
830 H(dd, dd) = 1;
831 invH(dd, dd) = 1;
832 }
833
834 a_res.resize(3);
835 int nb_gauss_pts = row_data.getN().size1();
836 commonData.valVel.resize(nb_gauss_pts);
837 commonData.jacVelRowPtr.resize(nb_gauss_pts);
838 commonData.jacVel.resize(nb_gauss_pts);
839
840 int nb_active_vars = 0;
841 for (int gg = 0; gg < nb_gauss_pts; gg++) {
842
843 if (gg == 0) {
844
845 trace_on(tAg);
846
847 for (int nn1 = 0; nn1 < 3; nn1++) { // 0
848 v[nn1] <<=
850 nb_active_vars++;
851 }
852 for (int nn1 = 0; nn1 < 3; nn1++) { // 3
853 dot_w[nn1] <<=
855 [gg][nn1];
856 nb_active_vars++;
857 }
859 .size() > 0) {
860 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3 = 6
861 for (int nn2 = 0; nn2 < 3; nn2++) {
862 h(nn1, nn2) <<=
864 nn1, nn2);
865 if (fieldDisp) {
866 if (nn1 == nn2) {
867 h(nn1, nn2) += 1;
868 }
869 }
870 nb_active_vars++;
871 }
872 }
873 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9
874 dot_W[nn1] <<=
876 .dataAtGaussPts["DOT_" + commonData.meshPositions][gg][nn1];
877 nb_active_vars++;
878 }
879 }
881 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9+3
882 for (int nn2 = 0; nn2 < 3; nn2++) {
883 H(nn1, nn2) <<=
885 nn2);
886 nb_active_vars++;
887 }
888 }
889 }
890 detH = 1;
891
895
899 auto t_invH = getFTensor2FromArray3by3(invH, FTensor::Number<0>(), 0);
900 auto t_dot_u =
901 FTensor::Tensor1<adouble *, 3>{&dot_u[0], &dot_u[1], &dot_u[2]};
902 auto t_dot_w =
903 FTensor::Tensor1<adouble *, 3>{&dot_w[0], &dot_w[1], &dot_w[2]};
904 auto t_dot_W =
905 FTensor::Tensor1<adouble *, 3>{&dot_W[0], &dot_W[1], &dot_W[2]};
906 auto t_v = FTensor::Tensor1<adouble *, 3>{&v[0], &v[1], &v[2]};
907 auto t_a_res =
908 FTensor::Tensor1<adouble *, 3>{&a_res[0], &a_res[1], &a_res[2]};
909
911 detH = determinantTensor3by3(H);
912 CHKERR invertTensor3by3(H, detH, invH);
913 t_F(i, j) = t_h(i, k) * t_invH(k, j);
914 } else {
915 t_F(i, j) = t_h(i, j);
916 }
917
918 t_dot_u(i) = t_dot_w(i) + t_F(i, j) * t_dot_W(j);
919 t_a_res(i) = t_v(i) - t_dot_u(i);
920 t_a_res(i) *= detH;
921
922 // dependant
923 VectorDouble &res = commonData.valVel[gg];
924 res.resize(3);
925 for (int rr = 0; rr < 3; rr++) {
926 a_res[rr] >>= res[rr];
927 }
928 trace_off();
929 }
930
931 active.resize(nb_active_vars);
932 int aa = 0;
933 for (int nn1 = 0; nn1 < 3; nn1++) {
934 active[aa++] =
936 }
937 for (int nn1 = 0; nn1 < 3; nn1++) {
938 active[aa++] =
940 .dataAtGaussPts["DOT_" + commonData.spatialPositions][gg][nn1];
941 }
943 0) {
944 for (int nn1 = 0; nn1 < 3; nn1++) {
945 for (int nn2 = 0; nn2 < 3; nn2++) {
946 if (fieldDisp && nn1 == nn2) {
947 active[aa++] =
949 nn1, nn2) +
950 1;
951 } else {
952 active[aa++] =
954 nn1, nn2);
955 }
956 }
957 }
958 for (int nn1 = 0; nn1 < 3; nn1++) {
959 active[aa++] =
961 .dataAtGaussPts["DOT_" + commonData.meshPositions][gg][nn1];
962 }
963 }
965 for (int nn1 = 0; nn1 < 3; nn1++) {
966 for (int nn2 = 0; nn2 < 3; nn2++) {
967 active[aa++] =
969 nn2);
970 }
971 }
972 }
973
974 if (!jAcobian) {
975 VectorDouble &res = commonData.valVel[gg];
976 if (gg > 0) {
977 res.resize(3);
978 int r;
979 r = ::function(tAg, 3, nb_active_vars, &active[0], &res[0]);
980 if (r != 3) {
981 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
982 "ADOL-C function evaluation with error");
983 }
984 }
985 double val = getVolume() * getGaussPts()(3, gg);
986 res *= val;
987 } else {
988 commonData.jacVelRowPtr[gg].resize(3);
989 commonData.jacVel[gg].resize(3, nb_active_vars);
990 for (int nn1 = 0; nn1 < 3; nn1++) {
991 (commonData.jacVelRowPtr[gg])[nn1] = &(commonData.jacVel[gg](nn1, 0));
992 }
993 int r;
994 r = jacobian(tAg, 3, nb_active_vars, &active[0],
995 &(commonData.jacVelRowPtr[gg])[0]);
996 if (r != 3) {
997 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
998 "ADOL-C function evaluation with error");
999 }
1000 double val = getVolume() * getGaussPts()(3, gg);
1001 commonData.jacVel[gg] *= val;
1002 // std::cerr << gg << " : " << commonData.jacVel[gg] << std::endl;
1003 }
1004 }
1005 }
1007}
1008
1010 const std::string field_name, BlockData &data, CommonData &common_data)
1013 dAta(data), commonData(common_data) {}
1014
1016 int row_side, EntityType row_type,
1017 EntitiesFieldData::EntData &row_data) {
1019
1020 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
1021 dAta.tEts.end()) {
1023 }
1024 int nb_dofs = row_data.getIndices().size();
1025 if (nb_dofs == 0)
1027
1028 auto base = row_data.getFTensor0N();
1029 int nb_base_functions = row_data.getN().size2();
1031
1032 {
1033
1034 nf.resize(nb_dofs);
1035 nf.clear();
1036
1037 for (unsigned int gg = 0; gg < row_data.getN().size1(); gg++) {
1038 FTensor::Tensor1<double *, 3> t_nf(&nf[0], &nf[1], &nf[2], 3);
1040 &commonData.valVel[gg][1],
1041 &commonData.valVel[gg][2]);
1042 int dd = 0;
1043 for (; dd < nb_dofs / 3; dd++) {
1044 t_nf(i) += base * res(i);
1045 ++base;
1046 ++t_nf;
1047 }
1048 for (; dd != nb_base_functions; dd++) {
1049 ++base;
1050 }
1051 }
1052
1053 if (row_data.getIndices().size() > 3 * row_data.getN().size2()) {
1054 SETERRQ(PETSC_COMM_SELF, 1, "data inconsistency");
1055 }
1056 CHKERR VecSetValues(getFEMethod()->ts_F, row_data.getIndices().size(),
1057 &row_data.getIndices()[0], &nf[0], ADD_VALUES);
1058 }
1060}
1061
1063 const std::string vel_field, const std::string field_name, BlockData &data,
1064 CommonData &common_data)
1065 : OpMassLhs_dM_dv(vel_field, field_name, data, common_data) {}
1066
1068 EntitiesFieldData::EntData &col_data, int gg) {
1070 int nb_col = col_data.getIndices().size();
1071 jac.clear();
1072 if (!nb_col)
1074 double *base_ptr = const_cast<double *>(&col_data.getN(gg)[0]);
1075 FTensor::Tensor0<double *> base(base_ptr, 1);
1076 FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
1077 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1078 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1080 &commonData.jacVel[gg](0, 0), &commonData.jacVel[gg](0, 1),
1081 &commonData.jacVel[gg](0, 2), &commonData.jacVel[gg](1, 0),
1082 &commonData.jacVel[gg](1, 1), &commonData.jacVel[gg](1, 2),
1083 &commonData.jacVel[gg](2, 0), &commonData.jacVel[gg](2, 1),
1084 &commonData.jacVel[gg](2, 2));
1087 for (int dd = 0; dd < nb_col / 3; dd++) {
1088 t_jac(i, j) += t_mass1(i, j) * base;
1089 ++base_ptr;
1090 ++t_jac;
1091 }
1092
1094}
1095
1097 const std::string vel_field, const std::string field_name, BlockData &data,
1098 CommonData &common_data)
1099 : OpVelocityLhs_dV_dv(vel_field, field_name, data, common_data) {}
1100
1102 EntitiesFieldData::EntData &col_data, int gg) {
1104 int nb_col = col_data.getIndices().size();
1105 jac.clear();
1106 if (!nb_col)
1108 double *base_ptr = const_cast<double *>(&col_data.getN(gg)[0]);
1109 FTensor::Tensor0<double *> base(base_ptr, 1);
1110 FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
1111 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1112 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1113 const int u = 3;
1115 &commonData.jacVel[gg](0, u + 0), &commonData.jacVel[gg](0, u + 1),
1116 &commonData.jacVel[gg](0, u + 2), &commonData.jacVel[gg](1, u + 0),
1117 &commonData.jacVel[gg](1, u + 1), &commonData.jacVel[gg](1, u + 2),
1118 &commonData.jacVel[gg](2, u + 0), &commonData.jacVel[gg](2, u + 1),
1119 &commonData.jacVel[gg](2, u + 2));
1123 if (commonData.dataAtGaussPts["DOT_" + commonData.meshPositions].size() ==
1124 0) {
1125
1126 for (int dd = 0; dd < nb_col / 3; dd++) {
1127 t_jac(i, j) += t_mass1(i, j) * base * getFEMethod()->ts_a;
1128 ++base_ptr;
1129 ++t_jac;
1130 }
1131 } else {
1132 double *diff_ptr =
1133 const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
1134 FTensor::Tensor1<double *, 3> diff(diff_ptr, &diff_ptr[1], &diff_ptr[2], 3);
1135 const int s = 3 + 3;
1137 // T* d000, T* d001, T* d002,
1138 // T* d010, T* d011, T* d012,
1139 // T* d020, T* d021, T* d022,
1140 // T* d100, T* d101, T* d102,
1141 // T* d110, T* d111, T* d112,
1142 // T* d120, T* d121, T* d122,
1143 // T* d200, T* d201, T* d202,
1144 // T* d210, T* d211, T* d212,
1145 // T* d220, T* d221, T* d222,
1146 &commonData.jacVel[gg](0, s + 0), &commonData.jacVel[gg](0, s + 1),
1147 &commonData.jacVel[gg](0, s + 2), &commonData.jacVel[gg](0, s + 3),
1148 &commonData.jacVel[gg](0, s + 4), &commonData.jacVel[gg](0, s + 5),
1149 &commonData.jacVel[gg](0, s + 6), &commonData.jacVel[gg](0, s + 7),
1150 &commonData.jacVel[gg](0, s + 8), &commonData.jacVel[gg](1, s + 0),
1151 &commonData.jacVel[gg](1, s + 1), &commonData.jacVel[gg](1, s + 2),
1152 &commonData.jacVel[gg](1, s + 3), &commonData.jacVel[gg](1, s + 4),
1153 &commonData.jacVel[gg](1, s + 5), &commonData.jacVel[gg](1, s + 6),
1154 &commonData.jacVel[gg](1, s + 7), &commonData.jacVel[gg](1, s + 8),
1155 &commonData.jacVel[gg](2, s + 0), &commonData.jacVel[gg](2, s + 1),
1156 &commonData.jacVel[gg](2, s + 2), &commonData.jacVel[gg](2, s + 3),
1157 &commonData.jacVel[gg](2, s + 4), &commonData.jacVel[gg](2, s + 5),
1158 &commonData.jacVel[gg](2, s + 6), &commonData.jacVel[gg](2, s + 7),
1159 &commonData.jacVel[gg](2, s + 8));
1160 for (int dd = 0; dd < nb_col / 3; dd++) {
1161 t_jac(i, j) += t_mass1(i, j) * base * getFEMethod()->ts_a;
1162 t_jac(i, j) += t_mass3(i, j, k) * diff(k);
1163 ++base_ptr;
1164 ++diff_ptr;
1165 ++t_jac;
1166 }
1167 }
1169}
1170
1172 const std::string vel_field, const std::string field_name, BlockData &data,
1173 CommonData &common_data)
1174 : OpVelocityLhs_dV_dv(vel_field, field_name, data, common_data) {}
1175
1177 EntitiesFieldData::EntData &col_data, int gg) {
1179 int nb_col = col_data.getIndices().size();
1180 jac.clear();
1181 if (!nb_col)
1183 double *base_ptr = const_cast<double *>(&col_data.getN(gg)[0]);
1184 FTensor::Tensor0<double *> base(base_ptr, 1);
1185 double *diff_ptr =
1186 const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
1187 FTensor::Tensor1<double *, 3> diff(diff_ptr, &diff_ptr[1], &diff_ptr[2], 3);
1188 FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
1189 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1190 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1191 const int u = 3 + 3 + 9;
1193 &commonData.jacVel[gg](0, u + 0), &commonData.jacVel[gg](0, u + 1),
1194 &commonData.jacVel[gg](0, u + 2), &commonData.jacVel[gg](1, u + 0),
1195 &commonData.jacVel[gg](1, u + 1), &commonData.jacVel[gg](1, u + 2),
1196 &commonData.jacVel[gg](2, u + 0), &commonData.jacVel[gg](2, u + 1),
1197 &commonData.jacVel[gg](2, u + 2));
1198 const int s = 3 + 3 + 9 + 3;
1200 // T* d000, T* d001, T* d002,
1201 // T* d010, T* d011, T* d012,
1202 // T* d020, T* d021, T* d022,
1203 // T* d100, T* d101, T* d102,
1204 // T* d110, T* d111, T* d112,
1205 // T* d120, T* d121, T* d122,
1206 // T* d200, T* d201, T* d202,
1207 // T* d210, T* d211, T* d212,
1208 // T* d220, T* d221, T* d222,
1209 &commonData.jacVel[gg](0, s + 0), &commonData.jacVel[gg](0, s + 1),
1210 &commonData.jacVel[gg](0, s + 2), &commonData.jacVel[gg](0, s + 3),
1211 &commonData.jacVel[gg](0, s + 4), &commonData.jacVel[gg](0, s + 5),
1212 &commonData.jacVel[gg](0, s + 6), &commonData.jacVel[gg](0, s + 7),
1213 &commonData.jacVel[gg](0, s + 8), &commonData.jacVel[gg](1, s + 0),
1214 &commonData.jacVel[gg](1, s + 1), &commonData.jacVel[gg](1, s + 2),
1215 &commonData.jacVel[gg](1, s + 3), &commonData.jacVel[gg](1, s + 4),
1216 &commonData.jacVel[gg](1, s + 5), &commonData.jacVel[gg](1, s + 6),
1217 &commonData.jacVel[gg](1, s + 7), &commonData.jacVel[gg](1, s + 8),
1218 &commonData.jacVel[gg](2, s + 0), &commonData.jacVel[gg](2, s + 1),
1219 &commonData.jacVel[gg](2, s + 2), &commonData.jacVel[gg](2, s + 3),
1220 &commonData.jacVel[gg](2, s + 4), &commonData.jacVel[gg](2, s + 5),
1221 &commonData.jacVel[gg](2, s + 6), &commonData.jacVel[gg](2, s + 7),
1222 &commonData.jacVel[gg](2, s + 8));
1226 for (int dd = 0; dd < nb_col / 3; dd++) {
1227 t_jac(i, j) += t_mass1(i, j) * base * getFEMethod()->ts_a;
1228 t_jac(i, j) += t_mass3(i, j, k) * diff(k);
1229 ++base_ptr;
1230 ++diff_ptr;
1231 ++t_jac;
1232 }
1234}
1235
1238 BlockData &data,
1239 CommonData &common_data, int tag,
1240 bool jacobian)
1243 dAta(data), commonData(common_data), tAg(tag), jAcobian(jacobian),
1244 fieldDisp(false) {}
1245
1248 int row_side, EntityType row_type,
1249 EntitiesFieldData::EntData &row_data) {
1251
1252 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
1253 dAta.tEts.end()) {
1255 }
1256
1257 // do it only once, no need to repeat this for edges,faces or tets
1258 if (row_type != MBVERTEX)
1260
1261 int nb_dofs = row_data.getIndices().size();
1262 if (nb_dofs == 0)
1264
1265 try {
1266
1267 a.resize(3);
1268 v.resize(3);
1269 g.resize(3, 3);
1270 G.resize(3, 3);
1271 h.resize(3, 3);
1272 F.resize(3, 3);
1273 H.resize(3, 3);
1274 H.clear();
1275 invH.resize(3, 3);
1276 invH.clear();
1277 for (int dd = 0; dd < 3; dd++) {
1278 H(dd, dd) = 1;
1279 invH(dd, dd) = 1;
1280 }
1281
1282 int nb_gauss_pts = row_data.getN().size1();
1283 commonData.valT.resize(nb_gauss_pts);
1284 commonData.jacTRowPtr.resize(nb_gauss_pts);
1285 commonData.jacT.resize(nb_gauss_pts);
1286
1287 int nb_active_vars = 0;
1288 for (int gg = 0; gg < nb_gauss_pts; gg++) {
1289
1290 if (gg == 0) {
1291
1292 trace_on(tAg);
1293
1294 for (int nn1 = 0; nn1 < 3; nn1++) { // 0
1295 a[nn1] <<=
1297 [gg][nn1];
1298 nb_active_vars++;
1299 }
1300
1301 for (int nn1 = 0; nn1 < 3; nn1++) { // 3
1302 v[nn1] <<=
1304 nb_active_vars++;
1305 }
1306 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3
1307 for (int nn2 = 0; nn2 < 3; nn2++) {
1308 g(nn1, nn2) <<=
1310 nn1, nn2);
1311 nb_active_vars++;
1312 }
1313 }
1314 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9
1315 for (int nn2 = 0; nn2 < 3; nn2++) {
1316 h(nn1, nn2) <<=
1318 nn2);
1319 nb_active_vars++;
1320 if (fieldDisp) {
1321 if (nn1 == nn2) {
1322 h(nn1, nn2) += 1;
1323 }
1324 }
1325 }
1326 }
1328 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9+9
1329 for (int nn2 = 0; nn2 < 3; nn2++) {
1330 H(nn1, nn2) <<=
1332 nn2);
1333 nb_active_vars++;
1334 }
1335 }
1336 }
1337 adouble detH;
1338 detH = 1;
1340 detH = determinantTensor3by3(H);
1341 }
1342 CHKERR invertTensor3by3(H, detH, invH);
1343
1347
1348 a_T.resize(3);
1349
1351 auto t_invH = getFTensor2FromArray3by3(invH, FTensor::Number<0>(), 0);
1352 auto t_F = getFTensor2FromArray3by3(F, FTensor::Number<0>(), 0);
1354 auto t_G = getFTensor2FromArray3by3(G, FTensor::Number<0>(), 0);
1355
1356 auto t_a = FTensor::Tensor1<adouble *, 3>{&a[0], &a[1], &a[2]};
1357 auto t_v = FTensor::Tensor1<adouble *, 3>{&v[0], &v[1], &v[2]};
1358 auto t_a_T = FTensor::Tensor1<adouble *, 3>{&a_T[0], &a_T[1], &a_T[2]};
1359
1360 t_F(i, j) = t_h(i, k) * t_invH(k, j);
1361 t_G(i, j) = t_g(i, k) * t_invH(k, j);
1362 t_a_T(i) = t_F(k, i) * t_a(k) + t_G(k, i) * t_v(k);
1363 const auto rho0 = dAta.rho0;
1364 t_a_T(i) = -rho0 * detH;
1365
1366 commonData.valT[gg].resize(3);
1367 for (int nn = 0; nn < 3; nn++) {
1368 a_T[nn] >>= (commonData.valT[gg])[nn];
1369 }
1370 trace_off();
1371 }
1372
1373 active.resize(nb_active_vars);
1374 int aa = 0;
1375 for (int nn1 = 0; nn1 < 3; nn1++) { // 0
1376 active[aa++] =
1378 .dataAtGaussPts["DOT_" + commonData.spatialVelocities][gg][nn1];
1379 }
1380
1381 for (int nn1 = 0; nn1 < 3; nn1++) { // 3
1382 active[aa++] =
1384 }
1385 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3
1386 for (int nn2 = 0; nn2 < 3; nn2++) {
1387 active[aa++] =
1389 nn2);
1390 }
1391 }
1392 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9
1393 for (int nn2 = 0; nn2 < 3; nn2++) {
1394 if (fieldDisp && nn1 == nn2) {
1395 active[aa++] =
1397 nn1, nn2) +
1398 1;
1399 } else {
1400 active[aa++] =
1402 nn2);
1403 }
1404 }
1405 }
1407 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9+9
1408 for (int nn2 = 0; nn2 < 3; nn2++) {
1409 active[aa++] =
1411 nn2);
1412 }
1413 }
1414 }
1415
1416 if (!jAcobian) {
1417 VectorDouble &res = commonData.valT[gg];
1418 if (gg > 0) {
1419 res.resize(3);
1420 int r;
1421 r = ::function(tAg, 3, nb_active_vars, &active[0], &res[0]);
1422 if (r != 3) { // function is locally analytic
1423 SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1424 "ADOL-C function evaluation with error r = %d", r);
1425 }
1426 }
1427 double val = getVolume() * getGaussPts()(3, gg);
1428 res *= val;
1429 } else {
1430 commonData.jacTRowPtr[gg].resize(3);
1431 commonData.jacT[gg].resize(3, nb_active_vars);
1432 for (int nn1 = 0; nn1 < 3; nn1++) {
1433 (commonData.jacTRowPtr[gg])[nn1] = &(commonData.jacT[gg](nn1, 0));
1434 }
1435 int r;
1436 r = jacobian(tAg, 3, nb_active_vars, &active[0],
1437 &(commonData.jacTRowPtr[gg])[0]);
1438 if (r != 3) {
1439 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1440 "ADOL-C function evaluation with error");
1441 }
1442 double val = getVolume() * getGaussPts()(3, gg);
1443 commonData.jacT[gg] *= val;
1444 }
1445 }
1446
1447 } catch (const std::exception &ex) {
1448 std::ostringstream ss;
1449 ss << "throw in method: " << ex.what() << std::endl;
1450 SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
1451 }
1452
1454}
1455
1458 BlockData &data,
1459 CommonData &common_data,
1460 Range *forcesonlyonentities_ptr)
1463 dAta(data), commonData(common_data) {
1464 if (forcesonlyonentities_ptr != NULL) {
1465 forcesOnlyOnEntities = *forcesonlyonentities_ptr;
1466 }
1467}
1468
1471 int row_side, EntityType row_type,
1472 EntitiesFieldData::EntData &row_data) {
1474
1475 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
1476 dAta.tEts.end()) {
1478 }
1479 int nb_dofs = row_data.getIndices().size();
1480 if (nb_dofs == 0)
1482
1483 try {
1484
1485 nf.resize(nb_dofs);
1486 nf.clear();
1487
1488 auto base = row_data.getFTensor0N();
1489 int nb_base_functions = row_data.getN().size2();
1491
1492 for (unsigned int gg = 0; gg < row_data.getN().size1(); gg++) {
1493 FTensor::Tensor1<double *, 3> t_nf(&nf[0], &nf[1], &nf[2], 3);
1495 &commonData.valT[gg][1],
1496 &commonData.valT[gg][2]);
1497 int dd = 0;
1498 for (; dd < nb_dofs / 3; dd++) {
1499 t_nf(i) += base * res(i);
1500 ++base;
1501 ++t_nf;
1502 }
1503 for (; dd != nb_base_functions; dd++) {
1504 ++base;
1505 }
1506 }
1507
1508 if (row_data.getIndices().size() > 3 * row_data.getN().size2()) {
1509 SETERRQ(PETSC_COMM_SELF, 1, "data inconsistency");
1510 }
1511 if (!forcesOnlyOnEntities.empty()) {
1512 VectorInt indices = row_data.getIndices();
1513 VectorDofs &dofs = row_data.getFieldDofs();
1514 VectorDofs::iterator dit = dofs.begin();
1515 for (int ii = 0; dit != dofs.end(); dit++, ii++) {
1516 if (forcesOnlyOnEntities.find((*dit)->getEnt()) ==
1517 forcesOnlyOnEntities.end()) {
1518 // std::cerr << **dit << std::endl;
1519 indices[ii] = -1;
1520 }
1521 }
1522 // std::cerr << indices << std::endl;
1523 CHKERR VecSetValues(getFEMethod()->ts_F, indices.size(), &indices[0],
1524 &nf[0], ADD_VALUES);
1525 } else {
1526 CHKERR VecSetValues(getFEMethod()->ts_F, row_data.getIndices().size(),
1527 &row_data.getIndices()[0], &nf[0], ADD_VALUES);
1528 }
1529
1530 } catch (const std::exception &ex) {
1531 std::ostringstream ss;
1532 ss << "throw in method: " << ex.what() << std::endl;
1533 SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
1534 }
1535
1537}
1538
1540 OpEshelbyDynamicMaterialMomentumLhs_dv(const std::string vel_field,
1541 const std::string field_name,
1542 BlockData &data,
1543 CommonData &common_data,
1544 Range *forcesonlyonentities_ptr)
1546 vel_field, field_name, data, common_data, forcesonlyonentities_ptr) {}
1547
1550 EntitiesFieldData::EntData &col_data, int gg) {
1552 int nb_col = col_data.getIndices().size();
1553 jac.clear();
1554 double *base_ptr = const_cast<double *>(&col_data.getN(gg)[0]);
1555 FTensor::Tensor0<double *> base(base_ptr, 1);
1556 double *diff_ptr =
1557 const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
1558 FTensor::Tensor1<double *, 3> diff(diff_ptr, &diff_ptr[1], &diff_ptr[2], 3);
1559 FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
1560 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1561 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1562 const int u = 3;
1564 &commonData.jacT[gg](0, u + 0), &commonData.jacT[gg](0, u + 1),
1565 &commonData.jacT[gg](0, u + 2), &commonData.jacT[gg](1, u + 0),
1566 &commonData.jacT[gg](1, u + 1), &commonData.jacT[gg](1, u + 2),
1567 &commonData.jacT[gg](2, u + 0), &commonData.jacT[gg](2, u + 1),
1568 &commonData.jacT[gg](2, u + 2));
1569 const int s = 3 + 3;
1571 // T* d000, T* d001, T* d002,
1572 // T* d010, T* d011, T* d012,
1573 // T* d020, T* d021, T* d022,
1574 // T* d100, T* d101, T* d102,
1575 // T* d110, T* d111, T* d112,
1576 // T* d120, T* d121, T* d122,
1577 // T* d200, T* d201, T* d202,
1578 // T* d210, T* d211, T* d212,
1579 // T* d220, T* d221, T* d222,
1580 &commonData.jacT[gg](0, s + 0), &commonData.jacT[gg](0, s + 1),
1581 &commonData.jacT[gg](0, s + 2), &commonData.jacT[gg](0, s + 3),
1582 &commonData.jacT[gg](0, s + 4), &commonData.jacT[gg](0, s + 5),
1583 &commonData.jacT[gg](0, s + 6), &commonData.jacT[gg](0, s + 7),
1584 &commonData.jacT[gg](0, s + 8), &commonData.jacT[gg](1, s + 0),
1585 &commonData.jacT[gg](1, s + 1), &commonData.jacT[gg](1, s + 2),
1586 &commonData.jacT[gg](1, s + 3), &commonData.jacT[gg](1, s + 4),
1587 &commonData.jacT[gg](1, s + 5), &commonData.jacT[gg](1, s + 6),
1588 &commonData.jacT[gg](1, s + 7), &commonData.jacT[gg](1, s + 8),
1589 &commonData.jacT[gg](2, s + 0), &commonData.jacT[gg](2, s + 1),
1590 &commonData.jacT[gg](2, s + 2), &commonData.jacT[gg](2, s + 3),
1591 &commonData.jacT[gg](2, s + 4), &commonData.jacT[gg](2, s + 5),
1592 &commonData.jacT[gg](2, s + 6), &commonData.jacT[gg](2, s + 7),
1593 &commonData.jacT[gg](2, s + 8));
1597 for (int dd = 0; dd < nb_col / 3; dd++) {
1598 t_jac(i, j) += t_mass1(i, j) * base * getFEMethod()->ts_a;
1599 t_jac(i, j) += t_mass3(i, j, k) * diff(k);
1600 ++base_ptr;
1601 ++diff_ptr;
1602 ++t_jac;
1603 }
1605}
1606
1608 OpEshelbyDynamicMaterialMomentumLhs_dx(const std::string vel_field,
1609 const std::string field_name,
1610 BlockData &data,
1611 CommonData &common_data,
1612 Range *forcesonlyonentities_ptr)
1614 vel_field, field_name, data, common_data, forcesonlyonentities_ptr) {}
1615
1618 EntitiesFieldData::EntData &col_data, int gg) {
1620 int nb_col = col_data.getIndices().size();
1621 jac.clear();
1622 double *diff_ptr =
1623 const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
1624 FTensor::Tensor1<double *, 3> diff(diff_ptr, &diff_ptr[1], &diff_ptr[2], 3);
1625 FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
1626 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1627 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1628 const int s = 3 + 3 + 9;
1630 // T* d000, T* d001, T* d002,
1631 // T* d010, T* d011, T* d012,
1632 // T* d020, T* d021, T* d022,
1633 // T* d100, T* d101, T* d102,
1634 // T* d110, T* d111, T* d112,
1635 // T* d120, T* d121, T* d122,
1636 // T* d200, T* d201, T* d202,
1637 // T* d210, T* d211, T* d212,
1638 // T* d220, T* d221, T* d222,
1639 &commonData.jacT[gg](0, s + 0), &commonData.jacT[gg](0, s + 1),
1640 &commonData.jacT[gg](0, s + 2), &commonData.jacT[gg](0, s + 3),
1641 &commonData.jacT[gg](0, s + 4), &commonData.jacT[gg](0, s + 5),
1642 &commonData.jacT[gg](0, s + 6), &commonData.jacT[gg](0, s + 7),
1643 &commonData.jacT[gg](0, s + 8), &commonData.jacT[gg](1, s + 0),
1644 &commonData.jacT[gg](1, s + 1), &commonData.jacT[gg](1, s + 2),
1645 &commonData.jacT[gg](1, s + 3), &commonData.jacT[gg](1, s + 4),
1646 &commonData.jacT[gg](1, s + 5), &commonData.jacT[gg](1, s + 6),
1647 &commonData.jacT[gg](1, s + 7), &commonData.jacT[gg](1, s + 8),
1648 &commonData.jacT[gg](2, s + 0), &commonData.jacT[gg](2, s + 1),
1649 &commonData.jacT[gg](2, s + 2), &commonData.jacT[gg](2, s + 3),
1650 &commonData.jacT[gg](2, s + 4), &commonData.jacT[gg](2, s + 5),
1651 &commonData.jacT[gg](2, s + 6), &commonData.jacT[gg](2, s + 7),
1652 &commonData.jacT[gg](2, s + 8));
1656 for (int dd = 0; dd < nb_col / 3; dd++) {
1657 t_jac(i, j) += t_mass3(i, j, k) * diff(k);
1658 ++diff_ptr;
1659 ++t_jac;
1660 }
1662}
1663
1665 OpEshelbyDynamicMaterialMomentumLhs_dX(const std::string vel_field,
1666 const std::string field_name,
1667 BlockData &data,
1668 CommonData &common_data,
1669 Range *forcesonlyonentities_ptr)
1671 vel_field, field_name, data, common_data, forcesonlyonentities_ptr) {}
1672
1675 EntitiesFieldData::EntData &col_data, int gg) {
1677 int nb_col = col_data.getIndices().size();
1678 jac.clear();
1679 double *diff_ptr =
1680 const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
1681 FTensor::Tensor1<double *, 3> diff(diff_ptr, &diff_ptr[1], &diff_ptr[2], 3);
1682 FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
1683 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1684 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1685 const int s = 3 + 3 + 9 + 9;
1687 // T* d000, T* d001, T* d002,
1688 // T* d010, T* d011, T* d012,
1689 // T* d020, T* d021, T* d022,
1690 // T* d100, T* d101, T* d102,
1691 // T* d110, T* d111, T* d112,
1692 // T* d120, T* d121, T* d122,
1693 // T* d200, T* d201, T* d202,
1694 // T* d210, T* d211, T* d212,
1695 // T* d220, T* d221, T* d222,
1696 &commonData.jacT[gg](0, s + 0), &commonData.jacT[gg](0, s + 1),
1697 &commonData.jacT[gg](0, s + 2), &commonData.jacT[gg](0, s + 3),
1698 &commonData.jacT[gg](0, s + 4), &commonData.jacT[gg](0, s + 5),
1699 &commonData.jacT[gg](0, s + 6), &commonData.jacT[gg](0, s + 7),
1700 &commonData.jacT[gg](0, s + 8), &commonData.jacT[gg](1, s + 0),
1701 &commonData.jacT[gg](1, s + 1), &commonData.jacT[gg](1, s + 2),
1702 &commonData.jacT[gg](1, s + 3), &commonData.jacT[gg](1, s + 4),
1703 &commonData.jacT[gg](1, s + 5), &commonData.jacT[gg](1, s + 6),
1704 &commonData.jacT[gg](1, s + 7), &commonData.jacT[gg](1, s + 8),
1705 &commonData.jacT[gg](2, s + 0), &commonData.jacT[gg](2, s + 1),
1706 &commonData.jacT[gg](2, s + 2), &commonData.jacT[gg](2, s + 3),
1707 &commonData.jacT[gg](2, s + 4), &commonData.jacT[gg](2, s + 5),
1708 &commonData.jacT[gg](2, s + 6), &commonData.jacT[gg](2, s + 7),
1709 &commonData.jacT[gg](2, s + 8));
1713 for (int dd = 0; dd < nb_col / 3; dd++) {
1714 t_jac(i, j) += t_mass3(i, j, k) * diff(k);
1715 ++diff_ptr;
1716 ++t_jac;
1717 }
1719}
1720
1722 MoFEM::Interface &m_field, TS _ts, const std::string velocity_field,
1723 const std::string spatial_position_field)
1724 : mField(m_field), tS(_ts), velocityField(velocity_field),
1725 spatialPositionField(spatial_position_field), jacobianLag(-1) {}
1726
1729
1730 switch (ts_ctx) {
1731 case CTX_TSSETIFUNCTION: {
1732 snes_f = ts_F;
1733 // FIXME: This global scattering because Kuu problem and Dynamic problem
1734 // not share partitions. Both problem should use the same partitioning to
1735 // resolve this problem.
1736 CHKERR mField.getInterface<VecManager>()->setGlobalGhostVector(
1737 problemPtr, COL, ts_u, INSERT_VALUES, SCATTER_REVERSE);
1738 CHKERR mField.getInterface<VecManager>()->setOtherGlobalGhostVector(
1739 problemPtr, velocityField, "DOT_" + velocityField, COL, ts_u_t,
1740 INSERT_VALUES, SCATTER_REVERSE);
1741 CHKERR mField.getInterface<VecManager>()->setOtherGlobalGhostVector(
1742 problemPtr, spatialPositionField, "DOT_" + spatialPositionField, COL,
1743 ts_u_t, INSERT_VALUES, SCATTER_REVERSE);
1744 break;
1745 }
1746 case CTX_TSSETIJACOBIAN: {
1747 snes_B = ts_B;
1748 break;
1749 }
1750 default:
1751 break;
1752 }
1753
1755}
1756
1759 //
1760 // SNES snes;
1761 // CHKERR TSGetSNES(tS,&snes);
1762 // CHKERR SNESSetLagJacobian(snes,jacobianLag);
1764}
1765
1768
1769 Range added_tets;
1771 mField, BLOCKSET | BODYFORCESSET, it)) {
1772 int id = it->getMeshsetId();
1773 EntityHandle meshset = it->getMeshset();
1774 CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET,
1775 setOfBlocks[id].tEts, true);
1776 added_tets.merge(setOfBlocks[id].tEts);
1777 Block_BodyForces mydata;
1778 CHKERR it->getAttributeDataStructure(mydata);
1779 setOfBlocks[id].rho0 = mydata.data.density;
1780 setOfBlocks[id].a0.resize(3);
1781 setOfBlocks[id].a0[0] = mydata.data.acceleration_x;
1782 setOfBlocks[id].a0[1] = mydata.data.acceleration_y;
1783 setOfBlocks[id].a0[2] = mydata.data.acceleration_z;
1784 // std::cerr << setOfBlocks[id].tEts << std::endl;
1785 }
1786
1788 mField, BLOCKSET | MAT_ELASTICSET, it)) {
1789 Mat_Elastic mydata;
1790 CHKERR it->getAttributeDataStructure(mydata);
1791 if (mydata.data.User1 == 0)
1792 continue;
1793 Range tets;
1794 EntityHandle meshset = it->getMeshset();
1795 CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets, true);
1796 tets = subtract(tets, added_tets);
1797 if (tets.empty())
1798 continue;
1799 int id = it->getMeshsetId();
1800 setOfBlocks[-id].tEts = tets;
1801 setOfBlocks[-id].rho0 = mydata.data.User1;
1802 setOfBlocks[-id].a0.resize(3);
1803 setOfBlocks[-id].a0[0] = mydata.data.User2;
1804 setOfBlocks[-id].a0[1] = mydata.data.User3;
1805 setOfBlocks[-id].a0[2] = mydata.data.User4;
1806 // std::cerr << setOfBlocks[id].tEts << std::endl;
1807 }
1808
1810}
1811
1813 MoFEM::Interface &m_field,
1814 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr) {
1816
1817 if (!block_sets_ptr)
1818 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1819 "Pointer to block of sets is null");
1820
1822 m_field, BLOCKSET | BODYFORCESSET, it)) {
1823 Block_BodyForces mydata;
1824 CHKERR it->getAttributeDataStructure(mydata);
1825 int id = it->getMeshsetId();
1826 auto &block_data = (*block_sets_ptr)[id];
1827 EntityHandle meshset = it->getMeshset();
1828 CHKERR m_field.get_moab().get_entities_by_dimension(meshset, 3,
1829 block_data.tEts, true);
1830 block_data.rho0 = mydata.data.density;
1831 block_data.a0.resize(3);
1832 block_data.a0[0] = mydata.data.acceleration_x;
1833 block_data.a0[1] = mydata.data.acceleration_y;
1834 block_data.a0[2] = mydata.data.acceleration_z;
1835 }
1836
1838}
1839
1841 string element_name, string velocity_field_name,
1842 string spatial_position_field_name, string material_position_field_name,
1843 bool ale, BitRefLevel bit) {
1845
1846 //
1847
1848 CHKERR mField.add_finite_element(element_name, MF_ZERO);
1850 velocity_field_name);
1852 velocity_field_name);
1854 velocity_field_name);
1856 element_name, spatial_position_field_name);
1858 element_name, spatial_position_field_name);
1860 element_name, spatial_position_field_name);
1861 if (mField.check_field(material_position_field_name)) {
1862 if (ale) {
1864 element_name, material_position_field_name);
1866 element_name, material_position_field_name);
1868 element_name, "DOT_" + material_position_field_name);
1869 }
1871 element_name, material_position_field_name);
1872 }
1874 element_name, "DOT_" + velocity_field_name);
1876 element_name, "DOT_" + spatial_position_field_name);
1877
1878 Range tets;
1879 if (bit.any()) {
1880 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1881 bit, BitRefLevel().set(), MBTET, tets);
1882 }
1883
1884 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
1885 for (; sit != setOfBlocks.end(); sit++) {
1886 Range add_tets = sit->second.tEts;
1887 if (!tets.empty()) {
1888 add_tets = intersect(add_tets, tets);
1889 }
1891 element_name);
1892 }
1893
1895}
1896
1898 string element_name, string velocity_field_name,
1899 string spatial_position_field_name, string material_position_field_name,
1900 bool ale, BitRefLevel bit) {
1902
1903 //
1904
1905 CHKERR mField.add_finite_element(element_name, MF_ZERO);
1907 velocity_field_name);
1909 velocity_field_name);
1911 velocity_field_name);
1913 element_name, spatial_position_field_name);
1915 element_name, spatial_position_field_name);
1916 if (mField.check_field(material_position_field_name)) {
1917 if (ale) {
1919 element_name, material_position_field_name);
1921 element_name, "DOT_" + material_position_field_name);
1922 }
1924 element_name, material_position_field_name);
1925 }
1927 element_name, "DOT_" + velocity_field_name);
1929 element_name, "DOT_" + spatial_position_field_name);
1930
1931 Range tets;
1932 if (bit.any()) {
1933 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1934 bit, BitRefLevel().set(), MBTET, tets);
1935 }
1936
1937 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
1938 for (; sit != setOfBlocks.end(); sit++) {
1939 Range add_tets = sit->second.tEts;
1940 if (!tets.empty()) {
1941 add_tets = intersect(add_tets, tets);
1942 }
1944 element_name);
1945 }
1946
1948}
1949
1951 string element_name, string velocity_field_name,
1952 string spatial_position_field_name, string material_position_field_name,
1953 bool ale, BitRefLevel bit, Range *intersected) {
1955
1956 //
1957
1958 CHKERR mField.add_finite_element(element_name, MF_ZERO);
1960 velocity_field_name);
1962 velocity_field_name);
1964 element_name, spatial_position_field_name);
1966 element_name, spatial_position_field_name);
1967 if (mField.check_field(material_position_field_name)) {
1968 if (ale) {
1970 element_name, material_position_field_name);
1972 element_name, material_position_field_name);
1974 element_name, "DOT_" + material_position_field_name);
1975 }
1977 element_name, material_position_field_name);
1978 }
1980 element_name, "DOT_" + velocity_field_name);
1982 element_name, "DOT_" + spatial_position_field_name);
1983
1984 Range tets;
1985 if (bit.any()) {
1986 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1987 bit, BitRefLevel().set(), MBTET, tets);
1988 }
1989 if (intersected != NULL) {
1990 if (tets.empty()) {
1991 tets = *intersected;
1992 } else {
1993 tets = intersect(*intersected, tets);
1994 }
1995 }
1996
1997 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
1998 for (; sit != setOfBlocks.end(); sit++) {
1999 Range add_tets = sit->second.tEts;
2000 if (!tets.empty()) {
2001 add_tets = intersect(add_tets, tets);
2002 }
2004 element_name);
2005 }
2006
2008}
2009
2011 string velocity_field_name, string spatial_position_field_name,
2012 string material_position_field_name, bool ale, bool linear) {
2014
2015 commonData.spatialPositions = spatial_position_field_name;
2016 commonData.meshPositions = material_position_field_name;
2017 commonData.spatialVelocities = velocity_field_name;
2018 commonData.lInear = linear;
2019
2020 // Rhs
2021 feMassRhs.getOpPtrVector().push_back(
2022 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2023 feMassRhs.getOpPtrVector().push_back(
2024 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2025 feMassRhs.getOpPtrVector().push_back(
2026 new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2027 feMassRhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2028 "DOT_" + spatial_position_field_name, commonData));
2029 if (mField.check_field(material_position_field_name)) {
2030 feMassRhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2031 material_position_field_name, commonData));
2032 if (ale) {
2033 feMassRhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2034 "DOT_" + material_position_field_name, commonData));
2035 } else {
2036 feMassRhs.meshPositionsFieldName = material_position_field_name;
2037 }
2038 }
2039 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
2040 for (; sit != setOfBlocks.end(); sit++) {
2041 feMassRhs.getOpPtrVector().push_back(
2042 new OpMassJacobian(spatial_position_field_name, sit->second, commonData,
2043 methodsOp, tAg, false));
2044 feMassRhs.getOpPtrVector().push_back(
2045 new OpMassRhs(spatial_position_field_name, sit->second, commonData));
2046 }
2047
2048 // Lhs
2049 feMassLhs.getOpPtrVector().push_back(
2050 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2051 feMassLhs.getOpPtrVector().push_back(
2052 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2053 feMassLhs.getOpPtrVector().push_back(
2054 new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2055 feMassLhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2056 "DOT_" + spatial_position_field_name, commonData));
2057 if (mField.check_field(material_position_field_name)) {
2058 feMassLhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2059 material_position_field_name, commonData));
2060 if (ale) {
2061 feMassLhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2062 "DOT_" + material_position_field_name, commonData));
2063 } else {
2064 feMassLhs.meshPositionsFieldName = material_position_field_name;
2065 }
2066 }
2067 sit = setOfBlocks.begin();
2068 for (; sit != setOfBlocks.end(); sit++) {
2069 feMassLhs.getOpPtrVector().push_back(
2070 new OpMassJacobian(spatial_position_field_name, sit->second, commonData,
2071 methodsOp, tAg, true));
2072 feMassLhs.getOpPtrVector().push_back(
2073 new OpMassLhs_dM_dv(spatial_position_field_name, velocity_field_name,
2074 sit->second, commonData));
2075 feMassLhs.getOpPtrVector().push_back(new OpMassLhs_dM_dx(
2076 spatial_position_field_name, spatial_position_field_name, sit->second,
2077 commonData));
2078 if (mField.check_field(material_position_field_name)) {
2079 if (ale) {
2080 feMassLhs.getOpPtrVector().push_back(new OpMassLhs_dM_dX(
2081 spatial_position_field_name, material_position_field_name,
2082 sit->second, commonData));
2083 } else {
2084 feMassLhs.meshPositionsFieldName = material_position_field_name;
2085 }
2086 }
2087 }
2088
2089 // Energy
2090 feEnergy.getOpPtrVector().push_back(
2091 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2092 feEnergy.getOpPtrVector().push_back(
2093 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2094 if (mField.check_field(material_position_field_name)) {
2095 feEnergy.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2096 material_position_field_name, commonData));
2097 feEnergy.meshPositionsFieldName = material_position_field_name;
2098 }
2099 sit = setOfBlocks.begin();
2100 for (; sit != setOfBlocks.end(); sit++) {
2101 feEnergy.getOpPtrVector().push_back(new OpEnergy(
2102 spatial_position_field_name, sit->second, commonData, feEnergy.V));
2103 }
2104
2106}
2107
2109 string velocity_field_name, string spatial_position_field_name,
2110 string material_position_field_name, bool ale) {
2112
2113 commonData.spatialPositions = spatial_position_field_name;
2114 commonData.meshPositions = material_position_field_name;
2115 commonData.spatialVelocities = velocity_field_name;
2116
2117 // Rhs
2118 feVelRhs.getOpPtrVector().push_back(
2119 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2120 feVelRhs.getOpPtrVector().push_back(
2121 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2122 feVelRhs.getOpPtrVector().push_back(
2123 new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2124 if (mField.check_field(material_position_field_name)) {
2125 feVelRhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2126 "DOT_" + spatial_position_field_name, commonData));
2127 feVelRhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2128 material_position_field_name, commonData));
2129 if (ale) {
2130 feVelRhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2131 "DOT_" + material_position_field_name, commonData));
2132 } else {
2133 feVelRhs.meshPositionsFieldName = material_position_field_name;
2134 }
2135 }
2136 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
2137 for (; sit != setOfBlocks.end(); sit++) {
2138 feVelRhs.getOpPtrVector().push_back(new OpVelocityJacobian(
2139 velocity_field_name, sit->second, commonData, tAg, false));
2140 feVelRhs.getOpPtrVector().push_back(
2141 new OpVelocityRhs(velocity_field_name, sit->second, commonData));
2142 }
2143
2144 // Lhs
2145 feVelLhs.getOpPtrVector().push_back(
2146 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2147 feVelLhs.getOpPtrVector().push_back(
2148 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2149 feVelLhs.getOpPtrVector().push_back(
2150 new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2151 if (mField.check_field(material_position_field_name)) {
2152 feVelLhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2153 "DOT_" + spatial_position_field_name, commonData));
2154 feVelLhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2155 material_position_field_name, commonData));
2156 if (ale) {
2157 feVelLhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2158 "DOT_" + material_position_field_name, commonData));
2159 } else {
2160 feVelLhs.meshPositionsFieldName = material_position_field_name;
2161 }
2162 }
2163 sit = setOfBlocks.begin();
2164 for (; sit != setOfBlocks.end(); sit++) {
2165 feVelLhs.getOpPtrVector().push_back(new OpVelocityJacobian(
2166 velocity_field_name, sit->second, commonData, tAg));
2167 feVelLhs.getOpPtrVector().push_back(new OpVelocityLhs_dV_dv(
2168 velocity_field_name, velocity_field_name, sit->second, commonData));
2169 feVelLhs.getOpPtrVector().push_back(new OpVelocityLhs_dV_dx(
2170 velocity_field_name, spatial_position_field_name, sit->second,
2171 commonData));
2172 if (mField.check_field(material_position_field_name)) {
2173 if (ale) {
2174 feVelLhs.getOpPtrVector().push_back(new OpVelocityLhs_dV_dX(
2175 velocity_field_name, material_position_field_name, sit->second,
2176 commonData));
2177 } else {
2178 feVelLhs.meshPositionsFieldName = material_position_field_name;
2179 }
2180 }
2181 }
2182
2184}
2185
2187 string velocity_field_name, string spatial_position_field_name,
2188 string material_position_field_name, Range *forces_on_entities_ptr) {
2190
2191 commonData.spatialPositions = spatial_position_field_name;
2192 commonData.meshPositions = material_position_field_name;
2193 commonData.spatialVelocities = velocity_field_name;
2194
2195 // Rhs
2196 feTRhs.getOpPtrVector().push_back(
2197 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2198 feTRhs.getOpPtrVector().push_back(
2199 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2200 feTRhs.getOpPtrVector().push_back(
2201 new OpGetCommonDataAtGaussPts(material_position_field_name, commonData));
2202 feTRhs.getOpPtrVector().push_back(
2203 new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2204
2205 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
2206 for (; sit != setOfBlocks.end(); sit++) {
2207 feTRhs.getOpPtrVector().push_back(
2209 material_position_field_name, sit->second, commonData, tAg, false));
2210 feTRhs.getOpPtrVector().push_back(new OpEshelbyDynamicMaterialMomentumRhs(
2211 material_position_field_name, sit->second, commonData,
2212 forces_on_entities_ptr));
2213 }
2214
2215 // Lhs
2216 feTLhs.getOpPtrVector().push_back(
2217 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2218 feTLhs.getOpPtrVector().push_back(
2219 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2220 feTLhs.getOpPtrVector().push_back(
2221 new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2222 if (mField.check_field(material_position_field_name)) {
2223 feTLhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2224 material_position_field_name, commonData));
2225 }
2226 sit = setOfBlocks.begin();
2227 for (; sit != setOfBlocks.end(); sit++) {
2228 feTLhs.getOpPtrVector().push_back(
2230 material_position_field_name, sit->second, commonData, tAg));
2231 feTLhs.getOpPtrVector().push_back(
2233 material_position_field_name, velocity_field_name, sit->second,
2234 commonData, forces_on_entities_ptr));
2235 feTLhs.getOpPtrVector().push_back(
2237 material_position_field_name, spatial_position_field_name,
2238 sit->second, commonData, forces_on_entities_ptr));
2239 feTLhs.getOpPtrVector().push_back(
2241 material_position_field_name, material_position_field_name,
2242 sit->second, commonData, forces_on_entities_ptr));
2243 }
2244
2246}
2247
2249 string velocity_field_name, string spatial_position_field_name,
2250 string material_position_field_name, bool linear) {
2252
2253 commonData.spatialPositions = spatial_position_field_name;
2254 commonData.meshPositions = material_position_field_name;
2255 commonData.spatialVelocities = velocity_field_name;
2256 commonData.lInear = linear;
2257
2258 // Rhs
2259 feMassRhs.getOpPtrVector().push_back(
2260 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2261 feMassRhs.getOpPtrVector().push_back(
2262 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2263 feMassRhs.getOpPtrVector().push_back(
2264 new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2265 if (mField.check_field(material_position_field_name)) {
2266 feMassRhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2267 material_position_field_name, commonData));
2268 feMassRhs.meshPositionsFieldName = material_position_field_name;
2269 }
2270 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
2271 for (; sit != setOfBlocks.end(); sit++) {
2272 feMassRhs.getOpPtrVector().push_back(
2273 new OpMassJacobian(spatial_position_field_name, sit->second, commonData,
2274 methodsOp, tAg, false));
2275 feMassRhs.getOpPtrVector().push_back(
2276 new OpMassRhs(spatial_position_field_name, sit->second, commonData));
2277 }
2278
2279 // Lhs
2280 feMassLhs.getOpPtrVector().push_back(
2281 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2282 feMassLhs.getOpPtrVector().push_back(
2283 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2284 feMassLhs.getOpPtrVector().push_back(
2285 new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2286 if (mField.check_field(material_position_field_name)) {
2287 feMassLhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2288 material_position_field_name, commonData));
2289 feMassLhs.meshPositionsFieldName = material_position_field_name;
2290 }
2291 sit = setOfBlocks.begin();
2292 for (; sit != setOfBlocks.end(); sit++) {
2293 feMassLhs.getOpPtrVector().push_back(
2294 new OpMassJacobian(spatial_position_field_name, sit->second, commonData,
2295 methodsOp, tAg, true));
2296 feMassLhs.getOpPtrVector().push_back(new OpMassLhs_dM_dv(
2297 spatial_position_field_name, spatial_position_field_name, sit->second,
2298 commonData));
2299 if (mField.check_field(material_position_field_name)) {
2300 feMassLhs.meshPositionsFieldName = material_position_field_name;
2301 }
2302 }
2303
2304 // Aux Lhs
2305 feMassAuxLhs.getOpPtrVector().push_back(
2306 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2307 feMassAuxLhs.getOpPtrVector().push_back(
2308 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2309 feMassAuxLhs.getOpPtrVector().push_back(
2310 new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2311 if (mField.check_field(material_position_field_name)) {
2312 feMassAuxLhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2313 material_position_field_name, commonData));
2314 feMassAuxLhs.meshPositionsFieldName = material_position_field_name;
2315 }
2316 sit = setOfBlocks.begin();
2317 for (; sit != setOfBlocks.end(); sit++) {
2318 feMassAuxLhs.getOpPtrVector().push_back(
2319 new OpMassJacobian(spatial_position_field_name, sit->second, commonData,
2320 methodsOp, tAg, true));
2321 feMassAuxLhs.getOpPtrVector().push_back(new OpMassLhs_dM_dx(
2322 spatial_position_field_name, spatial_position_field_name, sit->second,
2323 commonData));
2324 if (mField.check_field(material_position_field_name)) {
2325 feMassAuxLhs.meshPositionsFieldName = material_position_field_name;
2326 }
2327 }
2328
2329 // Energy E=0.5*rho*v*v
2330 feEnergy.getOpPtrVector().push_back(
2331 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2332 feEnergy.getOpPtrVector().push_back(
2333 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2334 if (mField.check_field(material_position_field_name)) {
2335 feEnergy.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2336 material_position_field_name, commonData));
2337 feEnergy.meshPositionsFieldName = material_position_field_name;
2338 }
2339 sit = setOfBlocks.begin();
2340 for (; sit != setOfBlocks.end(); sit++) {
2341 feEnergy.getOpPtrVector().push_back(new OpEnergy(
2342 spatial_position_field_name, sit->second, commonData, feEnergy.V));
2343 }
2344
2346}
2347
2350 if (iNitialized) {
2351
2352 CHKERR dEstroy();
2353 CHKERRABORT(PETSC_COMM_WORLD, ierr);
2354 }
2355}
2356
2359 if (!iNitialized) {
2360
2361#if PETSC_VERSION_GE(3, 5, 3)
2362 CHKERR MatCreateVecs(K, &u, &Ku);
2363 CHKERR MatCreateVecs(M, &v, &Mv);
2364#else
2365 CHKERR MatGetVecs(K, &u, &Ku);
2366 CHKERR MatGetVecs(M, &v, &Mv);
2367#endif
2368 CHKERR MatDuplicate(K, MAT_SHARE_NONZERO_PATTERN, &barK);
2369 iNitialized = true;
2370 }
2372}
2373
2376 if (iNitialized) {
2377
2378 CHKERR VecDestroy(&u);
2379 CHKERR VecDestroy(&Ku);
2380 CHKERR VecDestroy(&v);
2381 CHKERR VecDestroy(&Mv);
2382 CHKERR MatDestroy(&barK);
2383 iNitialized = false;
2384 }
2386}
2387
2390
2391 if (!initPC) {
2392 MPI_Comm comm;
2393 CHKERR PetscObjectGetComm((PetscObject)shellMat, &comm);
2394 CHKERR PCCreate(comm, &pC);
2395 initPC = true;
2396 }
2398}
2399
2402
2403 if (initPC) {
2404 CHKERR PCDestroy(&pC);
2405 initPC = false;
2406 }
2408}
2409
2411 MoFEM::Interface &m_field)
2412 : mField(m_field) {}
2413
2416
2417 if (ts_ctx != CTX_TSSETIFUNCTION) {
2418 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
2419 "It is used to residual of velocities");
2420 }
2421 if (!shellMatCtx->iNitialized) {
2422 CHKERR shellMatCtx->iNit();
2423 }
2424 // Note velocities calculate from displacements are stroed in shellMatCtx->u
2425 CHKERR VecScatterBegin(shellMatCtx->scatterU, ts_u_t, shellMatCtx->u,
2426 INSERT_VALUES, SCATTER_FORWARD);
2427 CHKERR VecScatterEnd(shellMatCtx->scatterU, ts_u_t, shellMatCtx->u,
2428 INSERT_VALUES, SCATTER_FORWARD);
2429 CHKERR VecScatterBegin(shellMatCtx->scatterV, ts_u, shellMatCtx->v,
2430 INSERT_VALUES, SCATTER_FORWARD);
2431 CHKERR VecScatterEnd(shellMatCtx->scatterV, ts_u, shellMatCtx->v,
2432 INSERT_VALUES, SCATTER_FORWARD);
2433 CHKERR VecAXPY(shellMatCtx->v, -1, shellMatCtx->u);
2434 CHKERR VecScatterBegin(shellMatCtx->scatterV, shellMatCtx->v, ts_F,
2435 ADD_VALUES, SCATTER_REVERSE);
2436 CHKERR VecScatterEnd(shellMatCtx->scatterV, shellMatCtx->v, ts_F, ADD_VALUES,
2437 SCATTER_REVERSE);
2438 // VecView(shellMatCtx->v,PETSC_VIEWER_STDOUT_WORLD);
2439
2441}
2442
2443
2444#ifdef __DIRICHLET_HPP__
2445
2446ConvectiveMassElement::ShellMatrixElement::ShellMatrixElement(
2447 MoFEM::Interface &m_field)
2448 : mField(m_field) {}
2449
2450MoFEMErrorCode ConvectiveMassElement::ShellMatrixElement::preProcess() {
2452
2453 if (ts_ctx != CTX_TSSETIJACOBIAN) {
2454 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
2455 "It is used to calculate shell matrix only");
2456 }
2457
2458 shellMatCtx->ts_a = ts_a;
2459 DirichletBcPtr->copyTs(*((TSMethod *)this)); // copy context for TSMethod
2460
2461 DirichletBcPtr->dIag = 1;
2462 DirichletBcPtr->ts_B = shellMatCtx->K;
2463 CHKERR MatZeroEntries(shellMatCtx->K);
2464 CHKERR mField.problem_basic_method_preProcess(problemName, *DirichletBcPtr);
2465 LoopsToDoType::iterator itk = loopK.begin();
2466 for (; itk != loopK.end(); itk++) {
2467 itk->second->copyTs(*((TSMethod *)this));
2468 itk->second->ts_B = shellMatCtx->K;
2469 CHKERR mField.loop_finite_elements(problemName, itk->first, *itk->second);
2470 }
2471 LoopsToDoType::iterator itam = loopAuxM.begin();
2472 for (; itam != loopAuxM.end(); itam++) {
2473 itam->second->copyTs(*((TSMethod *)this));
2474 itam->second->ts_B = shellMatCtx->K;
2475 CHKERR mField.loop_finite_elements(problemName, itam->first, *itam->second);
2476 }
2477 CHKERR mField.problem_basic_method_postProcess(problemName, *DirichletBcPtr);
2478 CHKERR MatAssemblyBegin(shellMatCtx->K, MAT_FINAL_ASSEMBLY);
2479 CHKERR MatAssemblyEnd(shellMatCtx->K, MAT_FINAL_ASSEMBLY);
2480
2481 DirichletBcPtr->dIag = 0;
2482 DirichletBcPtr->ts_B = shellMatCtx->M;
2483 CHKERR MatZeroEntries(shellMatCtx->M);
2484 // CHKERR mField.problem_basic_method_preProcess(problemName,*DirichletBcPtr);
2485 LoopsToDoType::iterator itm = loopM.begin();
2486 for (; itm != loopM.end(); itm++) {
2487 itm->second->copyTs(*((TSMethod *)this));
2488 itm->second->ts_B = shellMatCtx->M;
2489 CHKERR mField.loop_finite_elements(problemName, itm->first, *itm->second);
2490 }
2491 CHKERR mField.problem_basic_method_postProcess(problemName, *DirichletBcPtr);
2492 CHKERR MatAssemblyBegin(shellMatCtx->M, MAT_FINAL_ASSEMBLY);
2493 CHKERR MatAssemblyEnd(shellMatCtx->M, MAT_FINAL_ASSEMBLY);
2494
2495 // barK
2496 CHKERR MatZeroEntries(shellMatCtx->barK);
2497 CHKERR MatCopy(shellMatCtx->K, shellMatCtx->barK, SAME_NONZERO_PATTERN);
2498 CHKERR MatAXPY(shellMatCtx->barK, ts_a, shellMatCtx->M, SAME_NONZERO_PATTERN);
2499 CHKERR MatAssemblyBegin(shellMatCtx->barK, MAT_FINAL_ASSEMBLY);
2500 CHKERR MatAssemblyEnd(shellMatCtx->barK, MAT_FINAL_ASSEMBLY);
2501
2502 // Matrix View
2503 // MatView(shellMatCtx->barK,PETSC_VIEWER_DRAW_WORLD);//PETSC_VIEWER_STDOUT_WORLD);
2504 // std::string wait;
2505 // std::cin >> wait;
2506
2508}
2509
2510#endif //__DIRICHLET_HPP__
static Index< 'M', 3 > M
static Index< 'K', 3 > K
Operators and data structures for mass and convective mass element.
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr double a
@ COL
Definition: definitions.h:123
@ MF_ZERO
Definition: definitions.h:98
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ BODYFORCESSET
block name is "BODY_FORCES"
Definition: definitions.h:162
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
Definition: definitions.h:159
@ BLOCKSET
Definition: definitions.h:148
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
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 add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string &name_row)=0
set field col which finite element use
virtual bool check_field(const std::string &name) const =0
check if field is in database
virtual MoFEMErrorCode problem_basic_method_postProcess(const Problem *problem_ptr, BasicMethod &method, int verb=DEFAULT_VERBOSITY)=0
Set data for BasicMethod.
virtual MoFEMErrorCode loop_finite_elements(const std::string &problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
auto bit
set bit
constexpr double a0
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
UBlasVector< int > VectorInt
Definition: Types.hpp:67
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
auto createSmartVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
auto getFTensor2FromArray3by3(ublas::matrix< T, L, A > &data, const FTensor::Number< S > &, const size_t rr)
Definition: Templates.hpp:972
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
Definition: Templates.hpp:1204
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1193
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
const double r
rate factor
double A
double h
double rho
Definition: plastic.cpp:89
double H
Definition: plastic.cpp:91
constexpr auto field_name
constexpr double g
data for calculation inertia forces
common data used by volume elements
std::vector< std::vector< double * > > jacTRowPtr
std::map< std::string, std::vector< VectorDouble > > dataAtGaussPts
std::vector< std::vector< double * > > jacVelRowPtr
std::map< std::string, std::vector< MatrixDouble > > gradAtGaussPts
std::vector< std::vector< double * > > jacMassRowPtr
int getRule(int order)
it is used to calculate nb. of Gauss integration points
MyVolumeFE(MoFEM::Interface &m_field)
OpEnergy(const std::string field_name, BlockData &data, CommonData &common_data, SmartPetscObj< Vec > v)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpEshelbyDynamicMaterialMomentumJacobian(const std::string field_name, BlockData &data, CommonData &common_data, int tag, bool jacobian=true)
OpEshelbyDynamicMaterialMomentumLhs_dX(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data, Range *forcesonlyonentities_ptr)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
OpEshelbyDynamicMaterialMomentumLhs_dv(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data, Range *forcesonlyonentities_ptr)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
OpEshelbyDynamicMaterialMomentumLhs_dx(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data, Range *forcesonlyonentities_ptr)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpEshelbyDynamicMaterialMomentumRhs(const std::string field_name, BlockData &data, CommonData &common_data, Range *forcesonlyonentities_ptr)
OpGetCommonDataAtGaussPts(const std::string field_name, CommonData &common_data)
OpGetDataAtGaussPts(const std::string field_name, std::vector< VectorDouble > &values_at_gauss_pts, std::vector< MatrixDouble > &gardient_at_gauss_pts)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
operator calculating deformation gradient
OpMassJacobian(const std::string field_name, BlockData &data, CommonData &common_data, boost::ptr_vector< MethodForForceScaling > &methods_op, int tag, bool linear=false)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpMassLhs_dM_dX(const std::string field_name, const std::string col_field, BlockData &data, CommonData &common_data)
MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
OpMassLhs_dM_dv(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data, Range *forcesonlyonentities_ptr=NULL)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
OpMassLhs_dM_dx(const std::string field_name, const std::string col_field, BlockData &data, CommonData &common_data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpMassRhs(const std::string field_name, BlockData &data, CommonData &common_data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpVelocityJacobian(const std::string field_name, BlockData &data, CommonData &common_data, int tag, bool jacobian=true)
OpVelocityLhs_dV_dX(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
OpVelocityLhs_dV_dv(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
OpVelocityLhs_dV_dx(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
OpVelocityRhs(const std::string field_name, BlockData &data, CommonData &common_data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
MoFEMErrorCode preProcess()
Calculate inconsistency between approximation of velocities and velocities calculated from displaceme...
UpdateAndControl(MoFEM::Interface &m_field, TS _ts, const std::string velocity_field, const std::string spatial_position_field)
MoFEMErrorCode preProcess()
Scatter values from t_u_dt on the fields.
structure grouping operators and data used for calculation of mass (convective) element \ nonlinear_e...
ConvectiveMassElement(MoFEM::Interface &m_field, short int tag)
MoFEMErrorCode setVelocityOperators(string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool ale=false)
MyVolumeFE feEnergy
calculate kinetic energy
MoFEMErrorCode setShellMatrixMassOperators(string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool linear=false)
MyVolumeFE feVelRhs
calculate right hand side for tetrahedral elements
MoFEMErrorCode addVelocityElement(string element_name, string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool ale=false, BitRefLevel bit=BitRefLevel())
MoFEMErrorCode addConvectiveMassElement(string element_name, string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool ale=false, BitRefLevel bit=BitRefLevel())
MoFEMErrorCode setKinematicEshelbyOperators(string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", Range *forces_on_entities_ptr=NULL)
MoFEMErrorCode addEshelbyDynamicMaterialMomentum(string element_name, string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool ale=false, BitRefLevel bit=BitRefLevel(), Range *intersected=NULL)
MyVolumeFE feTRhs
calculate right hand side for tetrahedral elements
MyVolumeFE feMassRhs
calculate right hand side for tetrahedral elements
MyVolumeFE feTLhs
calculate left hand side for tetrahedral elements
MoFEMErrorCode setConvectiveMassOperators(string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool ale=false, bool linear=false)
boost::ptr_vector< MethodForForceScaling > methodsOp
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
MyVolumeFE feVelLhs
calculate left hand side for tetrahedral elements
static MoFEMErrorCode applyScale(const FEMethod *fe, boost::ptr_vector< MethodForForceScaling > &methods_op, VectorDouble &nf)
Managing BitRefLevels.
Body force data structure.
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode problem_basic_method_preProcess(const Problem *problem_ptr, BasicMethod &method, int verb=DEFAULT_VERBOSITY)=0
Set data for BasicMethod.
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
get dofs values
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1FieldData()
Return FTensor of rank 1, i.e. vector from filed data coeffinects.
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
const VectorInt & getIndices() const
Get global indices of dofs on entity.
structure to get information form mofem into EntitiesFieldData
Elastic material data structure.
data structure for TS (time stepping) context
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Vector manager is used to create vectors \mofem_vectors.
Definition: VecManager.hpp:23