v0.13.1
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
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
93 int side, EntityType type, EntitiesFieldData::EntData &data) {
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 //FIXME: there is error somewhere for nonlinear case
296 // test dam example with -is_linear 0
297 if (!lInear) {
298
299 t_F(i,j) = t_h(i,k)*t_invH(k,j);
300 t_a_res(i) *= rho0 * detH;
301 t_a_res(i) *= determinantTensor3by3(t_F);
302
303 } else {
304
305 t_a_res(i) *= rho0 * detH;
306
307 }
308
309 // dependant
311 res.resize(3);
312 for (int rr = 0; rr < 3; rr++) {
313 a_res[rr] >>= res[rr];
314 }
315
316 trace_off();
317 }
318
319 active.resize(nb_active_vars);
320 int aa = 0;
321 for (int nn1 = 0; nn1 < 3; nn1++) { // 0
322 active[aa++] = dot_spacial_vel[gg][nn1];
323 }
324 for (int nn1 = 0; nn1 < 3; nn1++) { // 3
325 for (int nn2 = 0; nn2 < 3; nn2++) {
326 if (fieldDisp && nn1 == nn2) {
327 active[aa++] = spatial_positions_grad[gg](nn1, nn2) + 1;
328 } else {
329 active[aa++] = spatial_positions_grad[gg](nn1, nn2);
330 }
331 }
332 }
334 0) {
335 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+9=12
336 for (int nn2 = 0; nn2 < 3; nn2++) {
337 active[aa++] = spatial_velocities_grad[gg](nn1, nn2);
338 }
339 }
340 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+9+9=21
341 active[aa++] = meshpos_vel[gg][nn1];
342 }
343 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+9+9+3=24
344 for (int nn2 = 0; nn2 < 3; nn2++) {
345 active[aa++] = mesh_positions_gradient[gg](nn1, nn2);
346 }
347 }
348 }
349
350 if (!jAcobian) {
352 if (gg > 0) {
353 res.resize(3);
354 int r;
355 r = ::function(tAg, 3, nb_active_vars, &active[0], &res[0]);
356 if (r != 3) { // function is locally analytic
357 SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
358 "ADOL-C function evaluation with error r = %d", r);
359 }
360 }
361 double val = getVolume() * getGaussPts()(3, gg);
362 res *= val;
363 // cout << "my res " << res << endl;
364 } else {
365 commonData.jacMassRowPtr[gg].resize(3);
366 commonData.jacMass[gg].resize(3, nb_active_vars);
367 for (int nn1 = 0; nn1 < 3; nn1++) {
368 (commonData.jacMassRowPtr[gg])[nn1] =
369 &(commonData.jacMass[gg](nn1, 0));
370 }
371 int r;
372 r = jacobian(tAg, 3, nb_active_vars, &active[0],
373 &(commonData.jacMassRowPtr[gg])[0]);
374 if (r != 3) {
375 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
376 "ADOL-C function evaluation with error");
377 }
378 double val = getVolume() * getGaussPts()(3, gg);
379 commonData.jacMass[gg] *= val;
380 }
381 }
382 }
383
385}
386
388 BlockData &data,
389 CommonData &common_data)
392 dAta(data), commonData(common_data) {}
393
395 int row_side, EntityType row_type,
396 EntitiesFieldData::EntData &row_data) {
398
399 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
400 dAta.tEts.end()) {
402 }
403 if (row_data.getIndices().size() == 0)
405 int nb_dofs = row_data.getIndices().size();
406
407 auto base = row_data.getFTensor0N();
408 int nb_base_functions = row_data.getN().size2();
409
410 {
411
412 nf.resize(nb_dofs);
413 nf.clear();
414
416
417 for (unsigned int gg = 0; gg < row_data.getN().size1(); gg++) {
418 FTensor::Tensor1<double *, 3> t_nf(&nf[0], &nf[1], &nf[2], 3);
420 &commonData.valMass[gg][1],
421 &commonData.valMass[gg][2]);
422 int dd = 0;
423 for (; dd < nb_dofs / 3; dd++) {
424 t_nf(i) += base * res(i);
425 ++base;
426 ++t_nf;
427 }
428 for (; dd != nb_base_functions; dd++) {
429 ++base;
430 }
431 }
432
433 if ((unsigned int)nb_dofs > 3 * row_data.getN().size2()) {
434 SETERRQ(PETSC_COMM_SELF, 1, "data inconsistency");
435 }
436 CHKERR VecSetValues(getFEMethod()->ts_F, nb_dofs, &row_data.getIndices()[0],
437 &nf[0], ADD_VALUES);
438 }
439
441}
442
444 const std::string vel_field, const std::string field_name, BlockData &data,
445 CommonData &common_data, Range *forcesonlyonentities_ptr)
447 vel_field, field_name,
449 dAta(data), commonData(common_data) {
450 sYmm = false;
451 if (forcesonlyonentities_ptr != NULL) {
452 forcesOnlyOnEntities = *forcesonlyonentities_ptr;
453 }
454}
455
457 EntitiesFieldData::EntData &col_data, int gg) {
459 int nb_col = col_data.getIndices().size();
460 jac.clear();
461 if (!nb_col)
466 FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
467 &jac(1, 0), &jac(1, 1), &jac(1, 2),
468 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
470 &commonData.jacMass[gg](0, 0), &commonData.jacMass[gg](0, 1),
471 &commonData.jacMass[gg](0, 2), &commonData.jacMass[gg](1, 0),
472 &commonData.jacMass[gg](1, 1), &commonData.jacMass[gg](1, 2),
473 &commonData.jacMass[gg](2, 0), &commonData.jacMass[gg](2, 1),
474 &commonData.jacMass[gg](2, 2));
475 double *base_ptr = const_cast<double *>(&col_data.getN(gg)[0]);
476 FTensor::Tensor0<double *> base(base_ptr, 1);
477 if (commonData.dataAtGaussPts["DOT_" + commonData.meshPositions].size() ==
478 0) {
479 for (int dd = 0; dd < nb_col / 3; dd++) {
480 t_jac(i, j) += t_mass1(i, j) * base * getFEMethod()->ts_a;
481 ++base;
482 ++t_jac;
483 }
484 } else {
485 const int s = 3 + 9;
487 // T* d000, T* d001, T* d002,
488 // T* d010, T* d011, T* d012,
489 // T* d020, T* d021, T* d022,
490 // T* d100, T* d101, T* d102,
491 // T* d110, T* d111, T* d112,
492 // T* d120, T* d121, T* d122,
493 // T* d200, T* d201, T* d202,
494 // T* d210, T* d211, T* d212,
495 // T* d220, T* d221, T* d222,
496 &commonData.jacMass[gg](0, s + 0), &commonData.jacMass[gg](0, s + 1),
497 &commonData.jacMass[gg](0, s + 2), &commonData.jacMass[gg](0, s + 3),
498 &commonData.jacMass[gg](0, s + 4), &commonData.jacMass[gg](0, s + 5),
499 &commonData.jacMass[gg](0, s + 6), &commonData.jacMass[gg](0, s + 7),
500 &commonData.jacMass[gg](0, s + 8), &commonData.jacMass[gg](1, s + 0),
501 &commonData.jacMass[gg](1, s + 1), &commonData.jacMass[gg](1, s + 2),
502 &commonData.jacMass[gg](1, s + 3), &commonData.jacMass[gg](1, s + 4),
503 &commonData.jacMass[gg](1, s + 5), &commonData.jacMass[gg](1, s + 6),
504 &commonData.jacMass[gg](1, s + 7), &commonData.jacMass[gg](1, s + 8),
505 &commonData.jacMass[gg](2, s + 0), &commonData.jacMass[gg](2, s + 1),
506 &commonData.jacMass[gg](2, s + 2), &commonData.jacMass[gg](2, s + 3),
507 &commonData.jacMass[gg](2, s + 4), &commonData.jacMass[gg](2, s + 5),
508 &commonData.jacMass[gg](2, s + 6), &commonData.jacMass[gg](2, s + 7),
509 &commonData.jacMass[gg](2, s + 8));
510
511 double *diff_ptr =
512 const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
513 FTensor::Tensor1<double *, 3> diff(diff_ptr, &diff_ptr[1], &diff_ptr[2], 3);
514 for (int dd = 0; dd < nb_col / 3; dd++) {
515 t_jac(i, j) += t_mass1(i, j) * base * getFEMethod()->ts_a;
516 t_jac(i, j) += t_mass3(i, j, k) * diff(k);
517 ++base;
518 ++diff;
519 ++t_jac;
520 }
521 }
523}
524
526 int row_side, int col_side, EntityType row_type, EntityType col_type,
528 EntitiesFieldData::EntData &col_data) {
530
531 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
532 dAta.tEts.end()) {
534 }
535
536 int nb_row = row_data.getIndices().size();
537 int nb_col = col_data.getIndices().size();
538 if (nb_row == 0)
540 if (nb_col == 0)
542
543 auto base = row_data.getFTensor0N();
544 int nb_base_functions = row_data.getN().size2();
545
546 {
547
548 k.resize(nb_row, nb_col);
549 k.clear();
550 jac.resize(3, nb_col);
551
552 for (unsigned int gg = 0; gg < row_data.getN().size1(); gg++) {
553
554 try {
555 CHKERR getJac(col_data, gg);
556 } catch (const std::exception &ex) {
557 std::ostringstream ss;
558 ss << "throw in method: " << ex.what() << std::endl;
559 SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
560 }
561
564
565 {
566 int dd1 = 0;
567 // integrate element stiffness matrix
568 for (; dd1 < nb_row / 3; dd1++) {
570 &jac(0, 0), &jac(0, 1), &jac(0, 2), &jac(1, 0), &jac(1, 1),
571 &jac(1, 2), &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
572 for (int dd2 = 0; dd2 < nb_col / 3; dd2++) {
574 &k(3 * dd1 + 0, 3 * dd2 + 0), &k(3 * dd1 + 0, 3 * dd2 + 1),
575 &k(3 * dd1 + 0, 3 * dd2 + 2), &k(3 * dd1 + 1, 3 * dd2 + 0),
576 &k(3 * dd1 + 1, 3 * dd2 + 1), &k(3 * dd1 + 1, 3 * dd2 + 2),
577 &k(3 * dd1 + 2, 3 * dd2 + 0), &k(3 * dd1 + 2, 3 * dd2 + 1),
578 &k(3 * dd1 + 2, 3 * dd2 + 2));
579 t_k(i, j) += base * t_jac(i, j);
580 ++t_jac;
581 }
582 ++base;
583 // for(int rr1 = 0;rr1<3;rr1++) {
584 // for(int dd2 = 0;dd2<nb_col;dd2++) {
585 // k(3*dd1+rr1,dd2) += row_data.getN()(gg,dd1)*jac(rr1,dd2);
586 // }
587 // }
588 }
589 for (; dd1 != nb_base_functions; dd1++) {
590 ++base;
591 }
592 }
593 }
594
595 if (!forcesOnlyOnEntities.empty()) {
596 VectorInt indices = row_data.getIndices();
597 VectorDofs &dofs = row_data.getFieldDofs();
598 VectorDofs::iterator dit = dofs.begin();
599 for (int ii = 0; dit != dofs.end(); dit++, ii++) {
600 if (forcesOnlyOnEntities.find((*dit)->getEnt()) ==
601 forcesOnlyOnEntities.end()) {
602 indices[ii] = -1;
603 }
604 }
605 CHKERR MatSetValues(getFEMethod()->ts_B, nb_row, &indices[0], nb_col,
606 &col_data.getIndices()[0], &k(0, 0), ADD_VALUES);
607 } else {
608 CHKERR MatSetValues(getFEMethod()->ts_B, nb_row,
609 &row_data.getIndices()[0], nb_col,
610 &col_data.getIndices()[0], &k(0, 0), ADD_VALUES);
611 }
612 }
614}
615
617 const std::string field_name, const std::string col_field, BlockData &data,
618 CommonData &common_data)
619 : OpMassLhs_dM_dv(field_name, col_field, data, common_data) {}
620
622 EntitiesFieldData::EntData &col_data, int gg) {
627 int nb_col = col_data.getIndices().size();
628 jac.clear();
629 FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
630 &jac(1, 0), &jac(1, 1), &jac(1, 2),
631 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
632 const int s = 3;
634 // T* d000, T* d001, T* d002,
635 // T* d010, T* d011, T* d012,
636 // T* d020, T* d021, T* d022,
637 // T* d100, T* d101, T* d102,
638 // T* d110, T* d111, T* d112,
639 // T* d120, T* d121, T* d122,
640 // T* d200, T* d201, T* d202,
641 // T* d210, T* d211, T* d212,
642 // T* d220, T* d221, T* d222,
643 &commonData.jacMass[gg](0, s + 0), &commonData.jacMass[gg](0, s + 1),
644 &commonData.jacMass[gg](0, s + 2), &commonData.jacMass[gg](0, s + 3),
645 &commonData.jacMass[gg](0, s + 4), &commonData.jacMass[gg](0, s + 5),
646 &commonData.jacMass[gg](0, s + 6), &commonData.jacMass[gg](0, s + 7),
647 &commonData.jacMass[gg](0, s + 8), &commonData.jacMass[gg](1, s + 0),
648 &commonData.jacMass[gg](1, s + 1), &commonData.jacMass[gg](1, s + 2),
649 &commonData.jacMass[gg](1, s + 3), &commonData.jacMass[gg](1, s + 4),
650 &commonData.jacMass[gg](1, s + 5), &commonData.jacMass[gg](1, s + 6),
651 &commonData.jacMass[gg](1, s + 7), &commonData.jacMass[gg](1, s + 8),
652 &commonData.jacMass[gg](2, s + 0), &commonData.jacMass[gg](2, s + 1),
653 &commonData.jacMass[gg](2, s + 2), &commonData.jacMass[gg](2, s + 3),
654 &commonData.jacMass[gg](2, s + 4), &commonData.jacMass[gg](2, s + 5),
655 &commonData.jacMass[gg](2, s + 6), &commonData.jacMass[gg](2, s + 7),
656 &commonData.jacMass[gg](2, s + 8));
657 double *diff_ptr =
658 const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
659 FTensor::Tensor1<double *, 3> diff(diff_ptr, &diff_ptr[1], &diff_ptr[2], 3);
660 for (int dd = 0; dd < nb_col / 3; dd++) {
661 t_jac(i, j) += t_mass3(i, j, k) * diff(k);
662 ++diff;
663 ++t_jac;
664 }
666}
667
669 const std::string field_name, const std::string col_field, BlockData &data,
670 CommonData &common_data)
671 : OpMassLhs_dM_dv(field_name, col_field, data, common_data) {}
672
674 EntitiesFieldData::EntData &col_data, int gg) {
676 int nb_col = col_data.getIndices().size();
677 jac.clear();
678 double *base_ptr = const_cast<double *>(&col_data.getN(gg)[0]);
679 FTensor::Tensor0<double *> base(base_ptr, 1);
680 double *diff_ptr =
681 const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
682 FTensor::Tensor1<double *, 3> diff(diff_ptr, &diff_ptr[1], &diff_ptr[2], 3);
683 FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
684 &jac(1, 0), &jac(1, 1), &jac(1, 2),
685 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
686 const int u = 3 + 9 + 9;
688 &commonData.jacMass[gg](0, u + 0), &commonData.jacMass[gg](0, u + 1),
689 &commonData.jacMass[gg](0, u + 2), &commonData.jacMass[gg](1, u + 0),
690 &commonData.jacMass[gg](1, u + 1), &commonData.jacMass[gg](1, u + 2),
691 &commonData.jacMass[gg](2, u + 0), &commonData.jacMass[gg](2, u + 1),
692 &commonData.jacMass[gg](2, u + 2));
693 const int s = 3 + 9 + 9 + 3;
695 // T* d000, T* d001, T* d002,
696 // T* d010, T* d011, T* d012,
697 // T* d020, T* d021, T* d022,
698 // T* d100, T* d101, T* d102,
699 // T* d110, T* d111, T* d112,
700 // T* d120, T* d121, T* d122,
701 // T* d200, T* d201, T* d202,
702 // T* d210, T* d211, T* d212,
703 // T* d220, T* d221, T* d222,
704 &commonData.jacMass[gg](0, s + 0), &commonData.jacMass[gg](0, s + 1),
705 &commonData.jacMass[gg](0, s + 2), &commonData.jacMass[gg](0, s + 3),
706 &commonData.jacMass[gg](0, s + 4), &commonData.jacMass[gg](0, s + 5),
707 &commonData.jacMass[gg](0, s + 6), &commonData.jacMass[gg](0, s + 7),
708 &commonData.jacMass[gg](0, s + 8), &commonData.jacMass[gg](1, s + 0),
709 &commonData.jacMass[gg](1, s + 1), &commonData.jacMass[gg](1, s + 2),
710 &commonData.jacMass[gg](1, s + 3), &commonData.jacMass[gg](1, s + 4),
711 &commonData.jacMass[gg](1, s + 5), &commonData.jacMass[gg](1, s + 6),
712 &commonData.jacMass[gg](1, s + 7), &commonData.jacMass[gg](1, s + 8),
713 &commonData.jacMass[gg](2, s + 0), &commonData.jacMass[gg](2, s + 1),
714 &commonData.jacMass[gg](2, s + 2), &commonData.jacMass[gg](2, s + 3),
715 &commonData.jacMass[gg](2, s + 4), &commonData.jacMass[gg](2, s + 5),
716 &commonData.jacMass[gg](2, s + 6), &commonData.jacMass[gg](2, s + 7),
717 &commonData.jacMass[gg](2, s + 8));
721 for (int dd = 0; dd < nb_col / 3; dd++) {
722 t_jac(i, j) += t_mass1(i, j) * base * getFEMethod()->ts_a;
723 t_jac(i, j) += t_mass3(i, j, k) * diff(k);
724 ++base_ptr;
725 ++diff_ptr;
726 ++t_jac;
727 }
729}
730
732 BlockData &data,
733 CommonData &common_data,
737 dAta(data), commonData(common_data), V(v, true),
738 lInear(commonData.lInear) {}
739
741 int row_side, EntityType row_type,
742 EntitiesFieldData::EntData &row_data) {
744
745 if (row_type != MBVERTEX) {
747 }
748 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
749 dAta.tEts.end()) {
751 }
752
753 {
754 double energy = 0;
755 for (unsigned int gg = 0; gg < row_data.getN().size1(); gg++) {
756 double val = getVolume() * getGaussPts()(3, gg);
757 double rho0 = dAta.rho0;
758 double rho;
759 if (lInear) {
760 rho = rho0;
761 } else {
762 h.resize(3, 3);
763 noalias(h) =
766 .size() > 0) {
767 H.resize(3, 3);
768 noalias(H) =
770 auto detH = determinantTensor3by3(H);
771 invH.resize(3, 3);
772 CHKERR invertTensor3by3(H, detH, invH);
773 F.resize(3, 3);
774 noalias(F) = prod(h, invH);
775 } else {
776 F.resize(3, 3);
777 noalias(F) = h;
778 }
779 double detF = determinantTensor3by3(F);
780 rho = detF * rho0;
781 }
782 v.resize(3);
784 energy += 0.5 * (rho * val) * inner_prod(v, v);
785 }
786 CHKERR VecSetValue(V, 0, energy, ADD_VALUES);
787 }
788
790}
791
793 const std::string field_name, BlockData &data, CommonData &common_data,
794 int tag, bool jacobian)
797 dAta(data), commonData(common_data), tAg(tag), jAcobian(jacobian),
798 fieldDisp(false) {}
799
801 int row_side, EntityType row_type,
802 EntitiesFieldData::EntData &row_data) {
804
805 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
806 dAta.tEts.end()) {
808 }
809
810 // do it only once, no need to repeat this for edges,faces or tets
811 if (row_type != MBVERTEX)
813
814 int nb_dofs = row_data.getIndices().size();
815 if (nb_dofs == 0)
817
818 {
819
820 v.resize(3);
821 dot_w.resize(3);
822 h.resize(3, 3);
823 h.clear();
824 F.resize(3, 3);
825 dot_W.resize(3);
826 dot_W.clear();
827 H.resize(3, 3);
828 H.clear();
829 invH.resize(3, 3);
830 invH.clear();
831 dot_u.resize(3);
832 for (int dd = 0; dd < 3; dd++) {
833 H(dd, dd) = 1;
834 invH(dd, dd) = 1;
835 }
836
837 a_res.resize(3);
838 int nb_gauss_pts = row_data.getN().size1();
839 commonData.valVel.resize(nb_gauss_pts);
840 commonData.jacVelRowPtr.resize(nb_gauss_pts);
841 commonData.jacVel.resize(nb_gauss_pts);
842
843 int nb_active_vars = 0;
844 for (int gg = 0; gg < nb_gauss_pts; gg++) {
845
846 if (gg == 0) {
847
848 trace_on(tAg);
849
850 for (int nn1 = 0; nn1 < 3; nn1++) { // 0
851 v[nn1] <<=
853 nb_active_vars++;
854 }
855 for (int nn1 = 0; nn1 < 3; nn1++) { // 3
856 dot_w[nn1] <<=
858 [gg][nn1];
859 nb_active_vars++;
860 }
862 .size() > 0) {
863 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3 = 6
864 for (int nn2 = 0; nn2 < 3; nn2++) {
865 h(nn1, nn2) <<=
867 nn1, nn2);
868 if (fieldDisp) {
869 if (nn1 == nn2) {
870 h(nn1, nn2) += 1;
871 }
872 }
873 nb_active_vars++;
874 }
875 }
876 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9
877 dot_W[nn1] <<=
879 .dataAtGaussPts["DOT_" + commonData.meshPositions][gg][nn1];
880 nb_active_vars++;
881 }
882 }
884 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9+3
885 for (int nn2 = 0; nn2 < 3; nn2++) {
886 H(nn1, nn2) <<=
888 nn2);
889 nb_active_vars++;
890 }
891 }
892 }
893 detH = 1;
894
898
902 auto t_invH = getFTensor2FromArray3by3(invH, FTensor::Number<0>(), 0);
903 auto t_dot_u =
904 FTensor::Tensor1<adouble *, 3>{&dot_u[0], &dot_u[1], &dot_u[2]};
905 auto t_dot_w =
906 FTensor::Tensor1<adouble *, 3>{&dot_w[0], &dot_w[1], &dot_w[2]};
907 auto t_dot_W =
908 FTensor::Tensor1<adouble *, 3>{&dot_W[0], &dot_W[1], &dot_W[2]};
909 auto t_v = FTensor::Tensor1<adouble *, 3>{&v[0], &v[1], &v[2]};
910 auto t_a_res =
911 FTensor::Tensor1<adouble *, 3>{&a_res[0], &a_res[1], &a_res[2]};
912
914 detH = determinantTensor3by3(H);
915 CHKERR invertTensor3by3(H, detH, invH);
916 t_F(i, j) = t_h(i, k) * t_invH(k, j);
917 } else {
918 t_F(i, j) = t_h(i, j);
919 }
920
921 t_dot_u(i) = t_dot_w(i) + t_F(i, j) * t_dot_W(j);
922 t_a_res(i) = t_v(i) - t_dot_u(i);
923 t_a_res(i) *= detH;
924
925 // dependant
926 VectorDouble &res = commonData.valVel[gg];
927 res.resize(3);
928 for (int rr = 0; rr < 3; rr++) {
929 a_res[rr] >>= res[rr];
930 }
931 trace_off();
932 }
933
934 active.resize(nb_active_vars);
935 int aa = 0;
936 for (int nn1 = 0; nn1 < 3; nn1++) {
937 active[aa++] =
939 }
940 for (int nn1 = 0; nn1 < 3; nn1++) {
941 active[aa++] =
943 .dataAtGaussPts["DOT_" + commonData.spatialPositions][gg][nn1];
944 }
946 0) {
947 for (int nn1 = 0; nn1 < 3; nn1++) {
948 for (int nn2 = 0; nn2 < 3; nn2++) {
949 if (fieldDisp && nn1 == nn2) {
950 active[aa++] =
952 nn1, nn2) +
953 1;
954 } else {
955 active[aa++] =
957 nn1, nn2);
958 }
959 }
960 }
961 for (int nn1 = 0; nn1 < 3; nn1++) {
962 active[aa++] =
964 .dataAtGaussPts["DOT_" + commonData.meshPositions][gg][nn1];
965 }
966 }
968 for (int nn1 = 0; nn1 < 3; nn1++) {
969 for (int nn2 = 0; nn2 < 3; nn2++) {
970 active[aa++] =
972 nn2);
973 }
974 }
975 }
976
977 if (!jAcobian) {
978 VectorDouble &res = commonData.valVel[gg];
979 if (gg > 0) {
980 res.resize(3);
981 int r;
982 r = ::function(tAg, 3, nb_active_vars, &active[0], &res[0]);
983 if (r != 3) {
984 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
985 "ADOL-C function evaluation with error");
986 }
987 }
988 double val = getVolume() * getGaussPts()(3, gg);
989 res *= val;
990 } else {
991 commonData.jacVelRowPtr[gg].resize(3);
992 commonData.jacVel[gg].resize(3, nb_active_vars);
993 for (int nn1 = 0; nn1 < 3; nn1++) {
994 (commonData.jacVelRowPtr[gg])[nn1] = &(commonData.jacVel[gg](nn1, 0));
995 }
996 int r;
997 r = jacobian(tAg, 3, nb_active_vars, &active[0],
998 &(commonData.jacVelRowPtr[gg])[0]);
999 if (r != 3) {
1000 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1001 "ADOL-C function evaluation with error");
1002 }
1003 double val = getVolume() * getGaussPts()(3, gg);
1004 commonData.jacVel[gg] *= val;
1005 // std::cerr << gg << " : " << commonData.jacVel[gg] << std::endl;
1006 }
1007 }
1008 }
1010}
1011
1013 const std::string field_name, BlockData &data, CommonData &common_data)
1016 dAta(data), commonData(common_data) {}
1017
1019 int row_side, EntityType row_type,
1020 EntitiesFieldData::EntData &row_data) {
1022
1023 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
1024 dAta.tEts.end()) {
1026 }
1027 int nb_dofs = row_data.getIndices().size();
1028 if (nb_dofs == 0)
1030
1031 auto base = row_data.getFTensor0N();
1032 int nb_base_functions = row_data.getN().size2();
1034
1035 {
1036
1037 nf.resize(nb_dofs);
1038 nf.clear();
1039
1040 for (unsigned int gg = 0; gg < row_data.getN().size1(); gg++) {
1041 FTensor::Tensor1<double *, 3> t_nf(&nf[0], &nf[1], &nf[2], 3);
1043 &commonData.valVel[gg][1],
1044 &commonData.valVel[gg][2]);
1045 int dd = 0;
1046 for (; dd < nb_dofs / 3; dd++) {
1047 t_nf(i) += base * res(i);
1048 ++base;
1049 ++t_nf;
1050 }
1051 for (; dd != nb_base_functions; dd++) {
1052 ++base;
1053 }
1054 }
1055
1056 if (row_data.getIndices().size() > 3 * row_data.getN().size2()) {
1057 SETERRQ(PETSC_COMM_SELF, 1, "data inconsistency");
1058 }
1059 CHKERR VecSetValues(getFEMethod()->ts_F, row_data.getIndices().size(),
1060 &row_data.getIndices()[0], &nf[0], ADD_VALUES);
1061 }
1063}
1064
1066 const std::string vel_field, const std::string field_name, BlockData &data,
1067 CommonData &common_data)
1068 : OpMassLhs_dM_dv(vel_field, field_name, data, common_data) {}
1069
1071 EntitiesFieldData::EntData &col_data, int gg) {
1073 int nb_col = col_data.getIndices().size();
1074 jac.clear();
1075 if (!nb_col)
1077 double *base_ptr = const_cast<double *>(&col_data.getN(gg)[0]);
1078 FTensor::Tensor0<double *> base(base_ptr, 1);
1079 FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
1080 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1081 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1083 &commonData.jacVel[gg](0, 0), &commonData.jacVel[gg](0, 1),
1084 &commonData.jacVel[gg](0, 2), &commonData.jacVel[gg](1, 0),
1085 &commonData.jacVel[gg](1, 1), &commonData.jacVel[gg](1, 2),
1086 &commonData.jacVel[gg](2, 0), &commonData.jacVel[gg](2, 1),
1087 &commonData.jacVel[gg](2, 2));
1090 for (int dd = 0; dd < nb_col / 3; dd++) {
1091 t_jac(i, j) += t_mass1(i, j) * base;
1092 ++base_ptr;
1093 ++t_jac;
1094 }
1095
1097}
1098
1100 const std::string vel_field, const std::string field_name, BlockData &data,
1101 CommonData &common_data)
1102 : OpVelocityLhs_dV_dv(vel_field, field_name, data, common_data) {}
1103
1105 EntitiesFieldData::EntData &col_data, int gg) {
1107 int nb_col = col_data.getIndices().size();
1108 jac.clear();
1109 if (!nb_col)
1111 double *base_ptr = const_cast<double *>(&col_data.getN(gg)[0]);
1112 FTensor::Tensor0<double *> base(base_ptr, 1);
1113 FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
1114 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1115 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1116 const int u = 3;
1118 &commonData.jacVel[gg](0, u + 0), &commonData.jacVel[gg](0, u + 1),
1119 &commonData.jacVel[gg](0, u + 2), &commonData.jacVel[gg](1, u + 0),
1120 &commonData.jacVel[gg](1, u + 1), &commonData.jacVel[gg](1, u + 2),
1121 &commonData.jacVel[gg](2, u + 0), &commonData.jacVel[gg](2, u + 1),
1122 &commonData.jacVel[gg](2, u + 2));
1126 if (commonData.dataAtGaussPts["DOT_" + commonData.meshPositions].size() ==
1127 0) {
1128
1129 for (int dd = 0; dd < nb_col / 3; dd++) {
1130 t_jac(i, j) += t_mass1(i, j) * base * getFEMethod()->ts_a;
1131 ++base_ptr;
1132 ++t_jac;
1133 }
1134 } else {
1135 double *diff_ptr =
1136 const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
1137 FTensor::Tensor1<double *, 3> diff(diff_ptr, &diff_ptr[1], &diff_ptr[2], 3);
1138 const int s = 3 + 3;
1140 // T* d000, T* d001, T* d002,
1141 // T* d010, T* d011, T* d012,
1142 // T* d020, T* d021, T* d022,
1143 // T* d100, T* d101, T* d102,
1144 // T* d110, T* d111, T* d112,
1145 // T* d120, T* d121, T* d122,
1146 // T* d200, T* d201, T* d202,
1147 // T* d210, T* d211, T* d212,
1148 // T* d220, T* d221, T* d222,
1149 &commonData.jacVel[gg](0, s + 0), &commonData.jacVel[gg](0, s + 1),
1150 &commonData.jacVel[gg](0, s + 2), &commonData.jacVel[gg](0, s + 3),
1151 &commonData.jacVel[gg](0, s + 4), &commonData.jacVel[gg](0, s + 5),
1152 &commonData.jacVel[gg](0, s + 6), &commonData.jacVel[gg](0, s + 7),
1153 &commonData.jacVel[gg](0, s + 8), &commonData.jacVel[gg](1, s + 0),
1154 &commonData.jacVel[gg](1, s + 1), &commonData.jacVel[gg](1, s + 2),
1155 &commonData.jacVel[gg](1, s + 3), &commonData.jacVel[gg](1, s + 4),
1156 &commonData.jacVel[gg](1, s + 5), &commonData.jacVel[gg](1, s + 6),
1157 &commonData.jacVel[gg](1, s + 7), &commonData.jacVel[gg](1, s + 8),
1158 &commonData.jacVel[gg](2, s + 0), &commonData.jacVel[gg](2, s + 1),
1159 &commonData.jacVel[gg](2, s + 2), &commonData.jacVel[gg](2, s + 3),
1160 &commonData.jacVel[gg](2, s + 4), &commonData.jacVel[gg](2, s + 5),
1161 &commonData.jacVel[gg](2, s + 6), &commonData.jacVel[gg](2, s + 7),
1162 &commonData.jacVel[gg](2, s + 8));
1163 for (int dd = 0; dd < nb_col / 3; dd++) {
1164 t_jac(i, j) += t_mass1(i, j) * base * getFEMethod()->ts_a;
1165 t_jac(i, j) += t_mass3(i, j, k) * diff(k);
1166 ++base_ptr;
1167 ++diff_ptr;
1168 ++t_jac;
1169 }
1170 }
1172}
1173
1175 const std::string vel_field, const std::string field_name, BlockData &data,
1176 CommonData &common_data)
1177 : OpVelocityLhs_dV_dv(vel_field, field_name, data, common_data) {}
1178
1180 EntitiesFieldData::EntData &col_data, int gg) {
1182 int nb_col = col_data.getIndices().size();
1183 jac.clear();
1184 if (!nb_col)
1186 double *base_ptr = const_cast<double *>(&col_data.getN(gg)[0]);
1187 FTensor::Tensor0<double *> base(base_ptr, 1);
1188 double *diff_ptr =
1189 const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
1190 FTensor::Tensor1<double *, 3> diff(diff_ptr, &diff_ptr[1], &diff_ptr[2], 3);
1191 FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
1192 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1193 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1194 const int u = 3 + 3 + 9;
1196 &commonData.jacVel[gg](0, u + 0), &commonData.jacVel[gg](0, u + 1),
1197 &commonData.jacVel[gg](0, u + 2), &commonData.jacVel[gg](1, u + 0),
1198 &commonData.jacVel[gg](1, u + 1), &commonData.jacVel[gg](1, u + 2),
1199 &commonData.jacVel[gg](2, u + 0), &commonData.jacVel[gg](2, u + 1),
1200 &commonData.jacVel[gg](2, u + 2));
1201 const int s = 3 + 3 + 9 + 3;
1203 // T* d000, T* d001, T* d002,
1204 // T* d010, T* d011, T* d012,
1205 // T* d020, T* d021, T* d022,
1206 // T* d100, T* d101, T* d102,
1207 // T* d110, T* d111, T* d112,
1208 // T* d120, T* d121, T* d122,
1209 // T* d200, T* d201, T* d202,
1210 // T* d210, T* d211, T* d212,
1211 // T* d220, T* d221, T* d222,
1212 &commonData.jacVel[gg](0, s + 0), &commonData.jacVel[gg](0, s + 1),
1213 &commonData.jacVel[gg](0, s + 2), &commonData.jacVel[gg](0, s + 3),
1214 &commonData.jacVel[gg](0, s + 4), &commonData.jacVel[gg](0, s + 5),
1215 &commonData.jacVel[gg](0, s + 6), &commonData.jacVel[gg](0, s + 7),
1216 &commonData.jacVel[gg](0, s + 8), &commonData.jacVel[gg](1, s + 0),
1217 &commonData.jacVel[gg](1, s + 1), &commonData.jacVel[gg](1, s + 2),
1218 &commonData.jacVel[gg](1, s + 3), &commonData.jacVel[gg](1, s + 4),
1219 &commonData.jacVel[gg](1, s + 5), &commonData.jacVel[gg](1, s + 6),
1220 &commonData.jacVel[gg](1, s + 7), &commonData.jacVel[gg](1, s + 8),
1221 &commonData.jacVel[gg](2, s + 0), &commonData.jacVel[gg](2, s + 1),
1222 &commonData.jacVel[gg](2, s + 2), &commonData.jacVel[gg](2, s + 3),
1223 &commonData.jacVel[gg](2, s + 4), &commonData.jacVel[gg](2, s + 5),
1224 &commonData.jacVel[gg](2, s + 6), &commonData.jacVel[gg](2, s + 7),
1225 &commonData.jacVel[gg](2, s + 8));
1229 for (int dd = 0; dd < nb_col / 3; dd++) {
1230 t_jac(i, j) += t_mass1(i, j) * base * getFEMethod()->ts_a;
1231 t_jac(i, j) += t_mass3(i, j, k) * diff(k);
1232 ++base_ptr;
1233 ++diff_ptr;
1234 ++t_jac;
1235 }
1237}
1238
1241 BlockData &data,
1242 CommonData &common_data, int tag,
1243 bool jacobian)
1246 dAta(data), commonData(common_data), tAg(tag), jAcobian(jacobian),
1247 fieldDisp(false) {}
1248
1251 int row_side, EntityType row_type,
1252 EntitiesFieldData::EntData &row_data) {
1254
1255 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
1256 dAta.tEts.end()) {
1258 }
1259
1260 // do it only once, no need to repeat this for edges,faces or tets
1261 if (row_type != MBVERTEX)
1263
1264 int nb_dofs = row_data.getIndices().size();
1265 if (nb_dofs == 0)
1267
1268 try {
1269
1270 a.resize(3);
1271 v.resize(3);
1272 g.resize(3, 3);
1273 G.resize(3, 3);
1274 h.resize(3, 3);
1275 F.resize(3, 3);
1276 H.resize(3, 3);
1277 H.clear();
1278 invH.resize(3, 3);
1279 invH.clear();
1280 for (int dd = 0; dd < 3; dd++) {
1281 H(dd, dd) = 1;
1282 invH(dd, dd) = 1;
1283 }
1284
1285 int nb_gauss_pts = row_data.getN().size1();
1286 commonData.valT.resize(nb_gauss_pts);
1287 commonData.jacTRowPtr.resize(nb_gauss_pts);
1288 commonData.jacT.resize(nb_gauss_pts);
1289
1290 int nb_active_vars = 0;
1291 for (int gg = 0; gg < nb_gauss_pts; gg++) {
1292
1293 if (gg == 0) {
1294
1295 trace_on(tAg);
1296
1297 for (int nn1 = 0; nn1 < 3; nn1++) { // 0
1298 a[nn1] <<=
1300 [gg][nn1];
1301 nb_active_vars++;
1302 }
1303
1304 for (int nn1 = 0; nn1 < 3; nn1++) { // 3
1305 v[nn1] <<=
1307 nb_active_vars++;
1308 }
1309 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3
1310 for (int nn2 = 0; nn2 < 3; nn2++) {
1311 g(nn1, nn2) <<=
1313 nn1, nn2);
1314 nb_active_vars++;
1315 }
1316 }
1317 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9
1318 for (int nn2 = 0; nn2 < 3; nn2++) {
1319 h(nn1, nn2) <<=
1321 nn2);
1322 nb_active_vars++;
1323 if (fieldDisp) {
1324 if (nn1 == nn2) {
1325 h(nn1, nn2) += 1;
1326 }
1327 }
1328 }
1329 }
1331 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9+9
1332 for (int nn2 = 0; nn2 < 3; nn2++) {
1333 H(nn1, nn2) <<=
1335 nn2);
1336 nb_active_vars++;
1337 }
1338 }
1339 }
1340 adouble detH;
1341 detH = 1;
1343 detH = determinantTensor3by3(H);
1344 }
1345 CHKERR invertTensor3by3(H, detH, invH);
1346
1350
1351 a_T.resize(3);
1352
1354 auto t_invH = getFTensor2FromArray3by3(invH, FTensor::Number<0>(), 0);
1355 auto t_F = getFTensor2FromArray3by3(F, FTensor::Number<0>(), 0);
1358
1359 auto t_a = FTensor::Tensor1<adouble *, 3>{&a[0], &a[1], &a[2]};
1360 auto t_v = FTensor::Tensor1<adouble *, 3>{&v[0], &v[1], &v[2]};
1361 auto t_a_T = FTensor::Tensor1<adouble *, 3>{&a_T[0], &a_T[1], &a_T[2]};
1362
1363 t_F(i, j) = t_h(i, k) * t_invH(k, j);
1364 t_G(i, j) = t_g(i, k) * t_invH(k, j);
1365 t_a_T(i) = t_F(k, i) * t_a(k) + t_G(k, i) * t_v(k);
1366 const auto rho0 = dAta.rho0;
1367 t_a_T(i) = -rho0 * detH;
1368
1369 commonData.valT[gg].resize(3);
1370 for (int nn = 0; nn < 3; nn++) {
1371 a_T[nn] >>= (commonData.valT[gg])[nn];
1372 }
1373 trace_off();
1374 }
1375
1376 active.resize(nb_active_vars);
1377 int aa = 0;
1378 for (int nn1 = 0; nn1 < 3; nn1++) { // 0
1379 active[aa++] =
1381 .dataAtGaussPts["DOT_" + commonData.spatialVelocities][gg][nn1];
1382 }
1383
1384 for (int nn1 = 0; nn1 < 3; nn1++) { // 3
1385 active[aa++] =
1387 }
1388 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3
1389 for (int nn2 = 0; nn2 < 3; nn2++) {
1390 active[aa++] =
1392 nn2);
1393 }
1394 }
1395 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9
1396 for (int nn2 = 0; nn2 < 3; nn2++) {
1397 if (fieldDisp && nn1 == nn2) {
1398 active[aa++] =
1400 nn1, nn2) +
1401 1;
1402 } else {
1403 active[aa++] =
1405 nn2);
1406 }
1407 }
1408 }
1410 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9+9
1411 for (int nn2 = 0; nn2 < 3; nn2++) {
1412 active[aa++] =
1414 nn2);
1415 }
1416 }
1417 }
1418
1419 if (!jAcobian) {
1420 VectorDouble &res = commonData.valT[gg];
1421 if (gg > 0) {
1422 res.resize(3);
1423 int r;
1424 r = ::function(tAg, 3, nb_active_vars, &active[0], &res[0]);
1425 if (r != 3) { // function is locally analytic
1426 SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1427 "ADOL-C function evaluation with error r = %d", r);
1428 }
1429 }
1430 double val = getVolume() * getGaussPts()(3, gg);
1431 res *= val;
1432 } else {
1433 commonData.jacTRowPtr[gg].resize(3);
1434 commonData.jacT[gg].resize(3, nb_active_vars);
1435 for (int nn1 = 0; nn1 < 3; nn1++) {
1436 (commonData.jacTRowPtr[gg])[nn1] = &(commonData.jacT[gg](nn1, 0));
1437 }
1438 int r;
1439 r = jacobian(tAg, 3, nb_active_vars, &active[0],
1440 &(commonData.jacTRowPtr[gg])[0]);
1441 if (r != 3) {
1442 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1443 "ADOL-C function evaluation with error");
1444 }
1445 double val = getVolume() * getGaussPts()(3, gg);
1446 commonData.jacT[gg] *= val;
1447 }
1448 }
1449
1450 } catch (const std::exception &ex) {
1451 std::ostringstream ss;
1452 ss << "throw in method: " << ex.what() << std::endl;
1453 SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
1454 }
1455
1457}
1458
1461 BlockData &data,
1462 CommonData &common_data,
1463 Range *forcesonlyonentities_ptr)
1466 dAta(data), commonData(common_data) {
1467 if (forcesonlyonentities_ptr != NULL) {
1468 forcesOnlyOnEntities = *forcesonlyonentities_ptr;
1469 }
1470}
1471
1474 int row_side, EntityType row_type,
1475 EntitiesFieldData::EntData &row_data) {
1477
1478 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
1479 dAta.tEts.end()) {
1481 }
1482 int nb_dofs = row_data.getIndices().size();
1483 if (nb_dofs == 0)
1485
1486 try {
1487
1488 nf.resize(nb_dofs);
1489 nf.clear();
1490
1491 auto base = row_data.getFTensor0N();
1492 int nb_base_functions = row_data.getN().size2();
1494
1495 for (unsigned int gg = 0; gg < row_data.getN().size1(); gg++) {
1496 FTensor::Tensor1<double *, 3> t_nf(&nf[0], &nf[1], &nf[2], 3);
1498 &commonData.valT[gg][1],
1499 &commonData.valT[gg][2]);
1500 int dd = 0;
1501 for (; dd < nb_dofs / 3; dd++) {
1502 t_nf(i) += base * res(i);
1503 ++base;
1504 ++t_nf;
1505 }
1506 for (; dd != nb_base_functions; dd++) {
1507 ++base;
1508 }
1509 }
1510
1511 if (row_data.getIndices().size() > 3 * row_data.getN().size2()) {
1512 SETERRQ(PETSC_COMM_SELF, 1, "data inconsistency");
1513 }
1514 if (!forcesOnlyOnEntities.empty()) {
1515 VectorInt indices = row_data.getIndices();
1516 VectorDofs &dofs = row_data.getFieldDofs();
1517 VectorDofs::iterator dit = dofs.begin();
1518 for (int ii = 0; dit != dofs.end(); dit++, ii++) {
1519 if (forcesOnlyOnEntities.find((*dit)->getEnt()) ==
1520 forcesOnlyOnEntities.end()) {
1521 // std::cerr << **dit << std::endl;
1522 indices[ii] = -1;
1523 }
1524 }
1525 // std::cerr << indices << std::endl;
1526 CHKERR VecSetValues(getFEMethod()->ts_F, indices.size(), &indices[0],
1527 &nf[0], ADD_VALUES);
1528 } else {
1529 CHKERR VecSetValues(getFEMethod()->ts_F, row_data.getIndices().size(),
1530 &row_data.getIndices()[0], &nf[0], ADD_VALUES);
1531 }
1532
1533 } catch (const std::exception &ex) {
1534 std::ostringstream ss;
1535 ss << "throw in method: " << ex.what() << std::endl;
1536 SETERRQ(PETSC_COMM_SELF, 1, ss.str().c_str());
1537 }
1538
1540}
1541
1543 OpEshelbyDynamicMaterialMomentumLhs_dv(const std::string vel_field,
1544 const std::string field_name,
1545 BlockData &data,
1546 CommonData &common_data,
1547 Range *forcesonlyonentities_ptr)
1549 vel_field, field_name, data, common_data, forcesonlyonentities_ptr) {}
1550
1553 EntitiesFieldData::EntData &col_data, int gg) {
1555 int nb_col = col_data.getIndices().size();
1556 jac.clear();
1557 double *base_ptr = const_cast<double *>(&col_data.getN(gg)[0]);
1558 FTensor::Tensor0<double *> base(base_ptr, 1);
1559 double *diff_ptr =
1560 const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
1561 FTensor::Tensor1<double *, 3> diff(diff_ptr, &diff_ptr[1], &diff_ptr[2], 3);
1562 FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
1563 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1564 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1565 const int u = 3;
1567 &commonData.jacT[gg](0, u + 0), &commonData.jacT[gg](0, u + 1),
1568 &commonData.jacT[gg](0, u + 2), &commonData.jacT[gg](1, u + 0),
1569 &commonData.jacT[gg](1, u + 1), &commonData.jacT[gg](1, u + 2),
1570 &commonData.jacT[gg](2, u + 0), &commonData.jacT[gg](2, u + 1),
1571 &commonData.jacT[gg](2, u + 2));
1572 const int s = 3 + 3;
1574 // T* d000, T* d001, T* d002,
1575 // T* d010, T* d011, T* d012,
1576 // T* d020, T* d021, T* d022,
1577 // T* d100, T* d101, T* d102,
1578 // T* d110, T* d111, T* d112,
1579 // T* d120, T* d121, T* d122,
1580 // T* d200, T* d201, T* d202,
1581 // T* d210, T* d211, T* d212,
1582 // T* d220, T* d221, T* d222,
1583 &commonData.jacT[gg](0, s + 0), &commonData.jacT[gg](0, s + 1),
1584 &commonData.jacT[gg](0, s + 2), &commonData.jacT[gg](0, s + 3),
1585 &commonData.jacT[gg](0, s + 4), &commonData.jacT[gg](0, s + 5),
1586 &commonData.jacT[gg](0, s + 6), &commonData.jacT[gg](0, s + 7),
1587 &commonData.jacT[gg](0, s + 8), &commonData.jacT[gg](1, s + 0),
1588 &commonData.jacT[gg](1, s + 1), &commonData.jacT[gg](1, s + 2),
1589 &commonData.jacT[gg](1, s + 3), &commonData.jacT[gg](1, s + 4),
1590 &commonData.jacT[gg](1, s + 5), &commonData.jacT[gg](1, s + 6),
1591 &commonData.jacT[gg](1, s + 7), &commonData.jacT[gg](1, s + 8),
1592 &commonData.jacT[gg](2, s + 0), &commonData.jacT[gg](2, s + 1),
1593 &commonData.jacT[gg](2, s + 2), &commonData.jacT[gg](2, s + 3),
1594 &commonData.jacT[gg](2, s + 4), &commonData.jacT[gg](2, s + 5),
1595 &commonData.jacT[gg](2, s + 6), &commonData.jacT[gg](2, s + 7),
1596 &commonData.jacT[gg](2, s + 8));
1600 for (int dd = 0; dd < nb_col / 3; dd++) {
1601 t_jac(i, j) += t_mass1(i, j) * base * getFEMethod()->ts_a;
1602 t_jac(i, j) += t_mass3(i, j, k) * diff(k);
1603 ++base_ptr;
1604 ++diff_ptr;
1605 ++t_jac;
1606 }
1608}
1609
1611 OpEshelbyDynamicMaterialMomentumLhs_dx(const std::string vel_field,
1612 const std::string field_name,
1613 BlockData &data,
1614 CommonData &common_data,
1615 Range *forcesonlyonentities_ptr)
1617 vel_field, field_name, data, common_data, forcesonlyonentities_ptr) {}
1618
1621 EntitiesFieldData::EntData &col_data, int gg) {
1623 int nb_col = col_data.getIndices().size();
1624 jac.clear();
1625 double *diff_ptr =
1626 const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
1627 FTensor::Tensor1<double *, 3> diff(diff_ptr, &diff_ptr[1], &diff_ptr[2], 3);
1628 FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
1629 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1630 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1631 const int s = 3 + 3 + 9;
1633 // T* d000, T* d001, T* d002,
1634 // T* d010, T* d011, T* d012,
1635 // T* d020, T* d021, T* d022,
1636 // T* d100, T* d101, T* d102,
1637 // T* d110, T* d111, T* d112,
1638 // T* d120, T* d121, T* d122,
1639 // T* d200, T* d201, T* d202,
1640 // T* d210, T* d211, T* d212,
1641 // T* d220, T* d221, T* d222,
1642 &commonData.jacT[gg](0, s + 0), &commonData.jacT[gg](0, s + 1),
1643 &commonData.jacT[gg](0, s + 2), &commonData.jacT[gg](0, s + 3),
1644 &commonData.jacT[gg](0, s + 4), &commonData.jacT[gg](0, s + 5),
1645 &commonData.jacT[gg](0, s + 6), &commonData.jacT[gg](0, s + 7),
1646 &commonData.jacT[gg](0, s + 8), &commonData.jacT[gg](1, s + 0),
1647 &commonData.jacT[gg](1, s + 1), &commonData.jacT[gg](1, s + 2),
1648 &commonData.jacT[gg](1, s + 3), &commonData.jacT[gg](1, s + 4),
1649 &commonData.jacT[gg](1, s + 5), &commonData.jacT[gg](1, s + 6),
1650 &commonData.jacT[gg](1, s + 7), &commonData.jacT[gg](1, s + 8),
1651 &commonData.jacT[gg](2, s + 0), &commonData.jacT[gg](2, s + 1),
1652 &commonData.jacT[gg](2, s + 2), &commonData.jacT[gg](2, s + 3),
1653 &commonData.jacT[gg](2, s + 4), &commonData.jacT[gg](2, s + 5),
1654 &commonData.jacT[gg](2, s + 6), &commonData.jacT[gg](2, s + 7),
1655 &commonData.jacT[gg](2, s + 8));
1659 for (int dd = 0; dd < nb_col / 3; dd++) {
1660 t_jac(i, j) += t_mass3(i, j, k) * diff(k);
1661 ++diff_ptr;
1662 ++t_jac;
1663 }
1665}
1666
1668 OpEshelbyDynamicMaterialMomentumLhs_dX(const std::string vel_field,
1669 const std::string field_name,
1670 BlockData &data,
1671 CommonData &common_data,
1672 Range *forcesonlyonentities_ptr)
1674 vel_field, field_name, data, common_data, forcesonlyonentities_ptr) {}
1675
1678 EntitiesFieldData::EntData &col_data, int gg) {
1680 int nb_col = col_data.getIndices().size();
1681 jac.clear();
1682 double *diff_ptr =
1683 const_cast<double *>(&(col_data.getDiffN(gg, nb_col / 3)(0, 0)));
1684 FTensor::Tensor1<double *, 3> diff(diff_ptr, &diff_ptr[1], &diff_ptr[2], 3);
1685 FTensor::Tensor2<double *, 3, 3> t_jac(&jac(0, 0), &jac(0, 1), &jac(0, 2),
1686 &jac(1, 0), &jac(1, 1), &jac(1, 2),
1687 &jac(2, 0), &jac(2, 1), &jac(2, 2), 3);
1688 const int s = 3 + 3 + 9 + 9;
1690 // T* d000, T* d001, T* d002,
1691 // T* d010, T* d011, T* d012,
1692 // T* d020, T* d021, T* d022,
1693 // T* d100, T* d101, T* d102,
1694 // T* d110, T* d111, T* d112,
1695 // T* d120, T* d121, T* d122,
1696 // T* d200, T* d201, T* d202,
1697 // T* d210, T* d211, T* d212,
1698 // T* d220, T* d221, T* d222,
1699 &commonData.jacT[gg](0, s + 0), &commonData.jacT[gg](0, s + 1),
1700 &commonData.jacT[gg](0, s + 2), &commonData.jacT[gg](0, s + 3),
1701 &commonData.jacT[gg](0, s + 4), &commonData.jacT[gg](0, s + 5),
1702 &commonData.jacT[gg](0, s + 6), &commonData.jacT[gg](0, s + 7),
1703 &commonData.jacT[gg](0, s + 8), &commonData.jacT[gg](1, s + 0),
1704 &commonData.jacT[gg](1, s + 1), &commonData.jacT[gg](1, s + 2),
1705 &commonData.jacT[gg](1, s + 3), &commonData.jacT[gg](1, s + 4),
1706 &commonData.jacT[gg](1, s + 5), &commonData.jacT[gg](1, s + 6),
1707 &commonData.jacT[gg](1, s + 7), &commonData.jacT[gg](1, s + 8),
1708 &commonData.jacT[gg](2, s + 0), &commonData.jacT[gg](2, s + 1),
1709 &commonData.jacT[gg](2, s + 2), &commonData.jacT[gg](2, s + 3),
1710 &commonData.jacT[gg](2, s + 4), &commonData.jacT[gg](2, s + 5),
1711 &commonData.jacT[gg](2, s + 6), &commonData.jacT[gg](2, s + 7),
1712 &commonData.jacT[gg](2, s + 8));
1716 for (int dd = 0; dd < nb_col / 3; dd++) {
1717 t_jac(i, j) += t_mass3(i, j, k) * diff(k);
1718 ++diff_ptr;
1719 ++t_jac;
1720 }
1722}
1723
1725 MoFEM::Interface &m_field, TS _ts, const std::string velocity_field,
1726 const std::string spatial_position_field)
1727 : mField(m_field), tS(_ts), velocityField(velocity_field),
1728 spatialPositionField(spatial_position_field), jacobianLag(-1) {}
1729
1732
1733 switch (ts_ctx) {
1734 case CTX_TSSETIFUNCTION: {
1735 snes_f = ts_F;
1736 // FIXME: This global scattering because Kuu problem and Dynamic problem
1737 // not share partitions. Both problem should use the same partitioning to
1738 // resolve this problem.
1739 CHKERR mField.getInterface<VecManager>()->setGlobalGhostVector(
1740 problemPtr, COL, ts_u, INSERT_VALUES, SCATTER_REVERSE);
1741 CHKERR mField.getInterface<VecManager>()->setOtherGlobalGhostVector(
1742 problemPtr, velocityField, "DOT_" + velocityField, COL, ts_u_t,
1743 INSERT_VALUES, SCATTER_REVERSE);
1744 CHKERR mField.getInterface<VecManager>()->setOtherGlobalGhostVector(
1745 problemPtr, spatialPositionField, "DOT_" + spatialPositionField, COL,
1746 ts_u_t, INSERT_VALUES, SCATTER_REVERSE);
1747 break;
1748 }
1749 case CTX_TSSETIJACOBIAN: {
1750 snes_B = ts_B;
1751 break;
1752 }
1753 default:
1754 break;
1755 }
1756
1758}
1759
1762 //
1763 // SNES snes;
1764 // CHKERR TSGetSNES(tS,&snes);
1765 // CHKERR SNESSetLagJacobian(snes,jacobianLag);
1767}
1768
1771
1772 Range added_tets;
1774 mField, BLOCKSET | BODYFORCESSET, it)) {
1775 int id = it->getMeshsetId();
1776 EntityHandle meshset = it->getMeshset();
1777 CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET,
1778 setOfBlocks[id].tEts, true);
1779 added_tets.merge(setOfBlocks[id].tEts);
1780 Block_BodyForces mydata;
1781 CHKERR it->getAttributeDataStructure(mydata);
1782 setOfBlocks[id].rho0 = mydata.data.density;
1783 setOfBlocks[id].a0.resize(3);
1784 setOfBlocks[id].a0[0] = mydata.data.acceleration_x;
1785 setOfBlocks[id].a0[1] = mydata.data.acceleration_y;
1786 setOfBlocks[id].a0[2] = mydata.data.acceleration_z;
1787 // std::cerr << setOfBlocks[id].tEts << std::endl;
1788 }
1789
1791 mField, BLOCKSET | MAT_ELASTICSET, it)) {
1792 Mat_Elastic mydata;
1793 CHKERR it->getAttributeDataStructure(mydata);
1794 if (mydata.data.User1 == 0)
1795 continue;
1796 Range tets;
1797 EntityHandle meshset = it->getMeshset();
1798 CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets, true);
1799 tets = subtract(tets, added_tets);
1800 if (tets.empty())
1801 continue;
1802 int id = it->getMeshsetId();
1803 setOfBlocks[-id].tEts = tets;
1804 setOfBlocks[-id].rho0 = mydata.data.User1;
1805 setOfBlocks[-id].a0.resize(3);
1806 setOfBlocks[-id].a0[0] = mydata.data.User2;
1807 setOfBlocks[-id].a0[1] = mydata.data.User3;
1808 setOfBlocks[-id].a0[2] = mydata.data.User4;
1809 // std::cerr << setOfBlocks[id].tEts << std::endl;
1810 }
1811
1813}
1814
1816 MoFEM::Interface &m_field,
1817 boost::shared_ptr<map<int, BlockData>> &block_sets_ptr) {
1819
1820 if (!block_sets_ptr)
1821 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1822 "Pointer to block of sets is null");
1823
1825 m_field, BLOCKSET | BODYFORCESSET, it)) {
1826 Block_BodyForces mydata;
1827 CHKERR it->getAttributeDataStructure(mydata);
1828 int id = it->getMeshsetId();
1829 auto &block_data = (*block_sets_ptr)[id];
1830 EntityHandle meshset = it->getMeshset();
1831 CHKERR m_field.get_moab().get_entities_by_dimension(meshset, 3,
1832 block_data.tEts, true);
1833 block_data.rho0 = mydata.data.density;
1834 block_data.a0.resize(3);
1835 block_data.a0[0] = mydata.data.acceleration_x;
1836 block_data.a0[1] = mydata.data.acceleration_y;
1837 block_data.a0[2] = mydata.data.acceleration_z;
1838 }
1839
1841}
1842
1844 string element_name, string velocity_field_name,
1845 string spatial_position_field_name, string material_position_field_name,
1846 bool ale, BitRefLevel bit) {
1848
1849 //
1850
1851 CHKERR mField.add_finite_element(element_name, MF_ZERO);
1853 velocity_field_name);
1855 velocity_field_name);
1857 velocity_field_name);
1859 element_name, spatial_position_field_name);
1861 element_name, spatial_position_field_name);
1863 element_name, spatial_position_field_name);
1864 if (mField.check_field(material_position_field_name)) {
1865 if (ale) {
1867 element_name, material_position_field_name);
1869 element_name, material_position_field_name);
1871 element_name, "DOT_" + material_position_field_name);
1872 }
1874 element_name, material_position_field_name);
1875 }
1877 element_name, "DOT_" + velocity_field_name);
1879 element_name, "DOT_" + spatial_position_field_name);
1880
1881 Range tets;
1882 if (bit.any()) {
1883 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1884 bit, BitRefLevel().set(), MBTET, tets);
1885 }
1886
1887 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
1888 for (; sit != setOfBlocks.end(); sit++) {
1889 Range add_tets = sit->second.tEts;
1890 if (!tets.empty()) {
1891 add_tets = intersect(add_tets, tets);
1892 }
1894 element_name);
1895 }
1896
1898}
1899
1901 string element_name, string velocity_field_name,
1902 string spatial_position_field_name, string material_position_field_name,
1903 bool ale, BitRefLevel bit) {
1905
1906 //
1907
1908 CHKERR mField.add_finite_element(element_name, MF_ZERO);
1910 velocity_field_name);
1912 velocity_field_name);
1914 velocity_field_name);
1916 element_name, spatial_position_field_name);
1918 element_name, spatial_position_field_name);
1919 if (mField.check_field(material_position_field_name)) {
1920 if (ale) {
1922 element_name, material_position_field_name);
1924 element_name, "DOT_" + material_position_field_name);
1925 }
1927 element_name, material_position_field_name);
1928 }
1930 element_name, "DOT_" + velocity_field_name);
1932 element_name, "DOT_" + spatial_position_field_name);
1933
1934 Range tets;
1935 if (bit.any()) {
1936 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1937 bit, BitRefLevel().set(), MBTET, tets);
1938 }
1939
1940 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
1941 for (; sit != setOfBlocks.end(); sit++) {
1942 Range add_tets = sit->second.tEts;
1943 if (!tets.empty()) {
1944 add_tets = intersect(add_tets, tets);
1945 }
1947 element_name);
1948 }
1949
1951}
1952
1954 string element_name, string velocity_field_name,
1955 string spatial_position_field_name, string material_position_field_name,
1956 bool ale, BitRefLevel bit, Range *intersected) {
1958
1959 //
1960
1961 CHKERR mField.add_finite_element(element_name, MF_ZERO);
1963 velocity_field_name);
1965 velocity_field_name);
1967 element_name, spatial_position_field_name);
1969 element_name, spatial_position_field_name);
1970 if (mField.check_field(material_position_field_name)) {
1971 if (ale) {
1973 element_name, material_position_field_name);
1975 element_name, material_position_field_name);
1977 element_name, "DOT_" + material_position_field_name);
1978 }
1980 element_name, material_position_field_name);
1981 }
1983 element_name, "DOT_" + velocity_field_name);
1985 element_name, "DOT_" + spatial_position_field_name);
1986
1987 Range tets;
1988 if (bit.any()) {
1989 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1990 bit, BitRefLevel().set(), MBTET, tets);
1991 }
1992 if (intersected != NULL) {
1993 if (tets.empty()) {
1994 tets = *intersected;
1995 } else {
1996 tets = intersect(*intersected, tets);
1997 }
1998 }
1999
2000 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
2001 for (; sit != setOfBlocks.end(); sit++) {
2002 Range add_tets = sit->second.tEts;
2003 if (!tets.empty()) {
2004 add_tets = intersect(add_tets, tets);
2005 }
2007 element_name);
2008 }
2009
2011}
2012
2014 string velocity_field_name, string spatial_position_field_name,
2015 string material_position_field_name, bool ale, bool linear) {
2017
2018 commonData.spatialPositions = spatial_position_field_name;
2019 commonData.meshPositions = material_position_field_name;
2020 commonData.spatialVelocities = velocity_field_name;
2021 commonData.lInear = linear;
2022
2023 // Rhs
2024 feMassRhs.getOpPtrVector().push_back(
2025 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2026 feMassRhs.getOpPtrVector().push_back(
2027 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2028 feMassRhs.getOpPtrVector().push_back(
2029 new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2030 feMassRhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2031 "DOT_" + spatial_position_field_name, commonData));
2032 if (mField.check_field(material_position_field_name)) {
2033 feMassRhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2034 material_position_field_name, commonData));
2035 if (ale) {
2036 feMassRhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2037 "DOT_" + material_position_field_name, commonData));
2038 } else {
2039 feMassRhs.meshPositionsFieldName = material_position_field_name;
2040 }
2041 }
2042 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
2043 for (; sit != setOfBlocks.end(); sit++) {
2044 feMassRhs.getOpPtrVector().push_back(
2045 new OpMassJacobian(spatial_position_field_name, sit->second, commonData,
2046 methodsOp, tAg, false));
2047 feMassRhs.getOpPtrVector().push_back(
2048 new OpMassRhs(spatial_position_field_name, sit->second, commonData));
2049 }
2050
2051 // Lhs
2052 feMassLhs.getOpPtrVector().push_back(
2053 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2054 feMassLhs.getOpPtrVector().push_back(
2055 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2056 feMassLhs.getOpPtrVector().push_back(
2057 new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2058 feMassLhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2059 "DOT_" + spatial_position_field_name, commonData));
2060 if (mField.check_field(material_position_field_name)) {
2061 feMassLhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2062 material_position_field_name, commonData));
2063 if (ale) {
2064 feMassLhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2065 "DOT_" + material_position_field_name, commonData));
2066 } else {
2067 feMassLhs.meshPositionsFieldName = material_position_field_name;
2068 }
2069 }
2070 sit = setOfBlocks.begin();
2071 for (; sit != setOfBlocks.end(); sit++) {
2072 feMassLhs.getOpPtrVector().push_back(
2073 new OpMassJacobian(spatial_position_field_name, sit->second, commonData,
2074 methodsOp, tAg, true));
2075 feMassLhs.getOpPtrVector().push_back(
2076 new OpMassLhs_dM_dv(spatial_position_field_name, velocity_field_name,
2077 sit->second, commonData));
2078 feMassLhs.getOpPtrVector().push_back(new OpMassLhs_dM_dx(
2079 spatial_position_field_name, spatial_position_field_name, sit->second,
2080 commonData));
2081 if (mField.check_field(material_position_field_name)) {
2082 if (ale) {
2083 feMassLhs.getOpPtrVector().push_back(new OpMassLhs_dM_dX(
2084 spatial_position_field_name, material_position_field_name,
2085 sit->second, commonData));
2086 } else {
2087 feMassLhs.meshPositionsFieldName = material_position_field_name;
2088 }
2089 }
2090 }
2091
2092 // Energy
2093 feEnergy.getOpPtrVector().push_back(
2094 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2095 feEnergy.getOpPtrVector().push_back(
2096 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2097 if (mField.check_field(material_position_field_name)) {
2098 feEnergy.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2099 material_position_field_name, commonData));
2100 feEnergy.meshPositionsFieldName = material_position_field_name;
2101 }
2102 sit = setOfBlocks.begin();
2103 for (; sit != setOfBlocks.end(); sit++) {
2104 feEnergy.getOpPtrVector().push_back(new OpEnergy(
2105 spatial_position_field_name, sit->second, commonData, feEnergy.V));
2106 }
2107
2109}
2110
2112 string velocity_field_name, string spatial_position_field_name,
2113 string material_position_field_name, bool ale) {
2115
2116 commonData.spatialPositions = spatial_position_field_name;
2117 commonData.meshPositions = material_position_field_name;
2118 commonData.spatialVelocities = velocity_field_name;
2119
2120 // Rhs
2121 feVelRhs.getOpPtrVector().push_back(
2122 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2123 feVelRhs.getOpPtrVector().push_back(
2124 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2125 feVelRhs.getOpPtrVector().push_back(
2126 new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2127 if (mField.check_field(material_position_field_name)) {
2128 feVelRhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2129 "DOT_" + spatial_position_field_name, commonData));
2130 feVelRhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2131 material_position_field_name, commonData));
2132 if (ale) {
2133 feVelRhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2134 "DOT_" + material_position_field_name, commonData));
2135 } else {
2136 feVelRhs.meshPositionsFieldName = material_position_field_name;
2137 }
2138 }
2139 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
2140 for (; sit != setOfBlocks.end(); sit++) {
2141 feVelRhs.getOpPtrVector().push_back(new OpVelocityJacobian(
2142 velocity_field_name, sit->second, commonData, tAg, false));
2143 feVelRhs.getOpPtrVector().push_back(
2144 new OpVelocityRhs(velocity_field_name, sit->second, commonData));
2145 }
2146
2147 // Lhs
2148 feVelLhs.getOpPtrVector().push_back(
2149 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2150 feVelLhs.getOpPtrVector().push_back(
2151 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2152 feVelLhs.getOpPtrVector().push_back(
2153 new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2154 if (mField.check_field(material_position_field_name)) {
2155 feVelLhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2156 "DOT_" + spatial_position_field_name, commonData));
2157 feVelLhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2158 material_position_field_name, commonData));
2159 if (ale) {
2160 feVelLhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2161 "DOT_" + material_position_field_name, commonData));
2162 } else {
2163 feVelLhs.meshPositionsFieldName = material_position_field_name;
2164 }
2165 }
2166 sit = setOfBlocks.begin();
2167 for (; sit != setOfBlocks.end(); sit++) {
2168 feVelLhs.getOpPtrVector().push_back(new OpVelocityJacobian(
2169 velocity_field_name, sit->second, commonData, tAg));
2170 feVelLhs.getOpPtrVector().push_back(new OpVelocityLhs_dV_dv(
2171 velocity_field_name, velocity_field_name, sit->second, commonData));
2172 feVelLhs.getOpPtrVector().push_back(new OpVelocityLhs_dV_dx(
2173 velocity_field_name, spatial_position_field_name, sit->second,
2174 commonData));
2175 if (mField.check_field(material_position_field_name)) {
2176 if (ale) {
2177 feVelLhs.getOpPtrVector().push_back(new OpVelocityLhs_dV_dX(
2178 velocity_field_name, material_position_field_name, sit->second,
2179 commonData));
2180 } else {
2181 feVelLhs.meshPositionsFieldName = material_position_field_name;
2182 }
2183 }
2184 }
2185
2187}
2188
2190 string velocity_field_name, string spatial_position_field_name,
2191 string material_position_field_name, Range *forces_on_entities_ptr) {
2193
2194 commonData.spatialPositions = spatial_position_field_name;
2195 commonData.meshPositions = material_position_field_name;
2196 commonData.spatialVelocities = velocity_field_name;
2197
2198 // Rhs
2199 feTRhs.getOpPtrVector().push_back(
2200 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2201 feTRhs.getOpPtrVector().push_back(
2202 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2203 feTRhs.getOpPtrVector().push_back(
2204 new OpGetCommonDataAtGaussPts(material_position_field_name, commonData));
2205 feTRhs.getOpPtrVector().push_back(
2206 new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2207
2208 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
2209 for (; sit != setOfBlocks.end(); sit++) {
2210 feTRhs.getOpPtrVector().push_back(
2212 material_position_field_name, sit->second, commonData, tAg, false));
2213 feTRhs.getOpPtrVector().push_back(new OpEshelbyDynamicMaterialMomentumRhs(
2214 material_position_field_name, sit->second, commonData,
2215 forces_on_entities_ptr));
2216 }
2217
2218 // Lhs
2219 feTLhs.getOpPtrVector().push_back(
2220 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2221 feTLhs.getOpPtrVector().push_back(
2222 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2223 feTLhs.getOpPtrVector().push_back(
2224 new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2225 if (mField.check_field(material_position_field_name)) {
2226 feTLhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2227 material_position_field_name, commonData));
2228 }
2229 sit = setOfBlocks.begin();
2230 for (; sit != setOfBlocks.end(); sit++) {
2231 feTLhs.getOpPtrVector().push_back(
2233 material_position_field_name, sit->second, commonData, tAg));
2234 feTLhs.getOpPtrVector().push_back(
2236 material_position_field_name, velocity_field_name, sit->second,
2237 commonData, forces_on_entities_ptr));
2238 feTLhs.getOpPtrVector().push_back(
2240 material_position_field_name, spatial_position_field_name,
2241 sit->second, commonData, forces_on_entities_ptr));
2242 feTLhs.getOpPtrVector().push_back(
2244 material_position_field_name, material_position_field_name,
2245 sit->second, commonData, forces_on_entities_ptr));
2246 }
2247
2249}
2250
2252 string velocity_field_name, string spatial_position_field_name,
2253 string material_position_field_name, bool linear) {
2255
2256 commonData.spatialPositions = spatial_position_field_name;
2257 commonData.meshPositions = material_position_field_name;
2258 commonData.spatialVelocities = velocity_field_name;
2259 commonData.lInear = linear;
2260
2261 // Rhs
2262 feMassRhs.getOpPtrVector().push_back(
2263 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2264 feMassRhs.getOpPtrVector().push_back(
2265 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2266 feMassRhs.getOpPtrVector().push_back(
2267 new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2268 if (mField.check_field(material_position_field_name)) {
2269 feMassRhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2270 material_position_field_name, commonData));
2271 feMassRhs.meshPositionsFieldName = material_position_field_name;
2272 }
2273 std::map<int, BlockData>::iterator sit = setOfBlocks.begin();
2274 for (; sit != setOfBlocks.end(); sit++) {
2275 feMassRhs.getOpPtrVector().push_back(
2276 new OpMassJacobian(spatial_position_field_name, sit->second, commonData,
2277 methodsOp, tAg, false));
2278 feMassRhs.getOpPtrVector().push_back(
2279 new OpMassRhs(spatial_position_field_name, sit->second, commonData));
2280 }
2281
2282 // Lhs
2283 feMassLhs.getOpPtrVector().push_back(
2284 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2285 feMassLhs.getOpPtrVector().push_back(
2286 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2287 feMassLhs.getOpPtrVector().push_back(
2288 new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2289 if (mField.check_field(material_position_field_name)) {
2290 feMassLhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2291 material_position_field_name, commonData));
2292 feMassLhs.meshPositionsFieldName = material_position_field_name;
2293 }
2294 sit = setOfBlocks.begin();
2295 for (; sit != setOfBlocks.end(); sit++) {
2296 feMassLhs.getOpPtrVector().push_back(
2297 new OpMassJacobian(spatial_position_field_name, sit->second, commonData,
2298 methodsOp, tAg, true));
2299 feMassLhs.getOpPtrVector().push_back(new OpMassLhs_dM_dv(
2300 spatial_position_field_name, spatial_position_field_name, sit->second,
2301 commonData));
2302 if (mField.check_field(material_position_field_name)) {
2303 feMassLhs.meshPositionsFieldName = material_position_field_name;
2304 }
2305 }
2306
2307 // Aux Lhs
2308 feMassAuxLhs.getOpPtrVector().push_back(
2309 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2310 feMassAuxLhs.getOpPtrVector().push_back(
2311 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2312 feMassAuxLhs.getOpPtrVector().push_back(
2313 new OpGetCommonDataAtGaussPts("DOT_" + velocity_field_name, commonData));
2314 if (mField.check_field(material_position_field_name)) {
2315 feMassAuxLhs.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2316 material_position_field_name, commonData));
2317 feMassAuxLhs.meshPositionsFieldName = material_position_field_name;
2318 }
2319 sit = setOfBlocks.begin();
2320 for (; sit != setOfBlocks.end(); sit++) {
2321 feMassAuxLhs.getOpPtrVector().push_back(
2322 new OpMassJacobian(spatial_position_field_name, sit->second, commonData,
2323 methodsOp, tAg, true));
2324 feMassAuxLhs.getOpPtrVector().push_back(new OpMassLhs_dM_dx(
2325 spatial_position_field_name, spatial_position_field_name, sit->second,
2326 commonData));
2327 if (mField.check_field(material_position_field_name)) {
2328 feMassAuxLhs.meshPositionsFieldName = material_position_field_name;
2329 }
2330 }
2331
2332 // Energy E=0.5*rho*v*v
2333 feEnergy.getOpPtrVector().push_back(
2334 new OpGetCommonDataAtGaussPts(velocity_field_name, commonData));
2335 feEnergy.getOpPtrVector().push_back(
2336 new OpGetCommonDataAtGaussPts(spatial_position_field_name, commonData));
2337 if (mField.check_field(material_position_field_name)) {
2338 feEnergy.getOpPtrVector().push_back(new OpGetCommonDataAtGaussPts(
2339 material_position_field_name, commonData));
2340 feEnergy.meshPositionsFieldName = material_position_field_name;
2341 }
2342 sit = setOfBlocks.begin();
2343 for (; sit != setOfBlocks.end(); sit++) {
2344 feEnergy.getOpPtrVector().push_back(new OpEnergy(
2345 spatial_position_field_name, sit->second, commonData, feEnergy.V));
2346 }
2347
2349}
2350
2353 if (iNitialized) {
2354
2355 CHKERR dEstroy();
2356 CHKERRABORT(PETSC_COMM_WORLD, ierr);
2357 }
2358}
2359
2362 if (!iNitialized) {
2363
2364#if PETSC_VERSION_GE(3, 5, 3)
2365 CHKERR MatCreateVecs(K, &u, &Ku);
2366 CHKERR MatCreateVecs(M, &v, &Mv);
2367#else
2368 CHKERR MatGetVecs(K, &u, &Ku);
2369 CHKERR MatGetVecs(M, &v, &Mv);
2370#endif
2371 CHKERR MatDuplicate(K, MAT_SHARE_NONZERO_PATTERN, &barK);
2372 iNitialized = true;
2373 }
2375}
2376
2379 if (iNitialized) {
2380
2381 CHKERR VecDestroy(&u);
2382 CHKERR VecDestroy(&Ku);
2383 CHKERR VecDestroy(&v);
2384 CHKERR VecDestroy(&Mv);
2385 CHKERR MatDestroy(&barK);
2386 iNitialized = false;
2387 }
2389}
2390
2393
2394 if (!initPC) {
2395 MPI_Comm comm;
2396 CHKERR PetscObjectGetComm((PetscObject)shellMat, &comm);
2397 CHKERR PCCreate(comm, &pC);
2398 initPC = true;
2399 }
2401}
2402
2405
2406 if (initPC) {
2407 CHKERR PCDestroy(&pC);
2408 initPC = false;
2409 }
2411}
2412
2414 MoFEM::Interface &m_field)
2415 : mField(m_field) {}
2416
2419
2420 if (ts_ctx != CTX_TSSETIFUNCTION) {
2421 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
2422 "It is used to residual of velocities");
2423 }
2424 if (!shellMatCtx->iNitialized) {
2425 CHKERR shellMatCtx->iNit();
2426 }
2427 // Note velocities calculate from displacements are stroed in shellMatCtx->u
2428 CHKERR VecScatterBegin(shellMatCtx->scatterU, ts_u_t, shellMatCtx->u,
2429 INSERT_VALUES, SCATTER_FORWARD);
2430 CHKERR VecScatterEnd(shellMatCtx->scatterU, ts_u_t, shellMatCtx->u,
2431 INSERT_VALUES, SCATTER_FORWARD);
2432 CHKERR VecScatterBegin(shellMatCtx->scatterV, ts_u, shellMatCtx->v,
2433 INSERT_VALUES, SCATTER_FORWARD);
2434 CHKERR VecScatterEnd(shellMatCtx->scatterV, ts_u, shellMatCtx->v,
2435 INSERT_VALUES, SCATTER_FORWARD);
2436 CHKERR VecAXPY(shellMatCtx->v, -1, shellMatCtx->u);
2437 CHKERR VecScatterBegin(shellMatCtx->scatterV, shellMatCtx->v, ts_F,
2438 ADD_VALUES, SCATTER_REVERSE);
2439 CHKERR VecScatterEnd(shellMatCtx->scatterV, shellMatCtx->v, ts_F, ADD_VALUES,
2440 SCATTER_REVERSE);
2441 // VecView(shellMatCtx->v,PETSC_VIEWER_STDOUT_WORLD);
2442
2444}
2445
2446
2447#ifdef __DIRICHLET_HPP__
2448
2449ConvectiveMassElement::ShellMatrixElement::ShellMatrixElement(
2450 MoFEM::Interface &m_field)
2451 : mField(m_field) {}
2452
2453MoFEMErrorCode ConvectiveMassElement::ShellMatrixElement::preProcess() {
2455
2456 if (ts_ctx != CTX_TSSETIJACOBIAN) {
2457 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
2458 "It is used to calculate shell matrix only");
2459 }
2460
2461 shellMatCtx->ts_a = ts_a;
2462 DirichletBcPtr->copyTs(*((TSMethod *)this)); // copy context for TSMethod
2463
2464 DirichletBcPtr->dIag = 1;
2465 DirichletBcPtr->ts_B = shellMatCtx->K;
2466 CHKERR MatZeroEntries(shellMatCtx->K);
2467 CHKERR mField.problem_basic_method_preProcess(problemName, *DirichletBcPtr);
2468 LoopsToDoType::iterator itk = loopK.begin();
2469 for (; itk != loopK.end(); itk++) {
2470 itk->second->copyTs(*((TSMethod *)this));
2471 itk->second->ts_B = shellMatCtx->K;
2472 CHKERR mField.loop_finite_elements(problemName, itk->first, *itk->second);
2473 }
2474 LoopsToDoType::iterator itam = loopAuxM.begin();
2475 for (; itam != loopAuxM.end(); itam++) {
2476 itam->second->copyTs(*((TSMethod *)this));
2477 itam->second->ts_B = shellMatCtx->K;
2478 CHKERR mField.loop_finite_elements(problemName, itam->first, *itam->second);
2479 }
2480 CHKERR mField.problem_basic_method_postProcess(problemName, *DirichletBcPtr);
2481 CHKERR MatAssemblyBegin(shellMatCtx->K, MAT_FINAL_ASSEMBLY);
2482 CHKERR MatAssemblyEnd(shellMatCtx->K, MAT_FINAL_ASSEMBLY);
2483
2484 DirichletBcPtr->dIag = 0;
2485 DirichletBcPtr->ts_B = shellMatCtx->M;
2486 CHKERR MatZeroEntries(shellMatCtx->M);
2487 // CHKERR mField.problem_basic_method_preProcess(problemName,*DirichletBcPtr);
2488 LoopsToDoType::iterator itm = loopM.begin();
2489 for (; itm != loopM.end(); itm++) {
2490 itm->second->copyTs(*((TSMethod *)this));
2491 itm->second->ts_B = shellMatCtx->M;
2492 CHKERR mField.loop_finite_elements(problemName, itm->first, *itm->second);
2493 }
2494 CHKERR mField.problem_basic_method_postProcess(problemName, *DirichletBcPtr);
2495 CHKERR MatAssemblyBegin(shellMatCtx->M, MAT_FINAL_ASSEMBLY);
2496 CHKERR MatAssemblyEnd(shellMatCtx->M, MAT_FINAL_ASSEMBLY);
2497
2498 // barK
2499 CHKERR MatZeroEntries(shellMatCtx->barK);
2500 CHKERR MatCopy(shellMatCtx->K, shellMatCtx->barK, SAME_NONZERO_PATTERN);
2501 CHKERR MatAXPY(shellMatCtx->barK, ts_a, shellMatCtx->M, SAME_NONZERO_PATTERN);
2502 CHKERR MatAssemblyBegin(shellMatCtx->barK, MAT_FINAL_ASSEMBLY);
2503 CHKERR MatAssemblyEnd(shellMatCtx->barK, MAT_FINAL_ASSEMBLY);
2504
2505 // Matrix View
2506 // MatView(shellMatCtx->barK,PETSC_VIEWER_DRAW_WORLD);//PETSC_VIEWER_STDOUT_WORLD);
2507 // std::string wait;
2508 // std::cin >> wait;
2509
2511}
2512
2513#endif //__DIRICHLET_HPP__
static Index< 'M', 3 > M
Operators and data structures for mass and convective mass element.
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr double a
@ COL
Definition: definitions.h:123
@ MF_ZERO
Definition: definitions.h:98
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ BODYFORCESSET
block name is "BODY_FORCES"
Definition: definitions.h:162
@ MAT_ELASTICSET
block name is "MAT_ELASTIC"
Definition: definitions.h:159
@ BLOCKSET
Definition: definitions.h:148
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
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
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
UBlasVector< int > VectorInt
Definition: Types.hpp:67
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
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:1091
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:1323
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1312
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
double h
double rho
Definition: plastic.cpp:98
double H
Definition: plastic.cpp:100
constexpr AssemblyType A
Definition: plastic.cpp:35
constexpr IntegrationType G
Definition: plastic.cpp:36
constexpr auto field_name
constexpr double g
data for calculation inertia forces
common data used by volume elements
std::vector< std::vector< double * > > jacTRowPtr
std::map< std::string, std::vector< VectorDouble > > dataAtGaussPts
std::vector< std::vector< double * > > jacVelRowPtr
std::map< std::string, std::vector< MatrixDouble > > gradAtGaussPts
std::vector< std::vector< double * > > jacMassRowPtr
int getRule(int order)
it is used to calculate nb. of Gauss integration points
MyVolumeFE(MoFEM::Interface &m_field)
OpEnergy(const std::string field_name, BlockData &data, CommonData &common_data, SmartPetscObj< Vec > v)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpEshelbyDynamicMaterialMomentumJacobian(const std::string field_name, BlockData &data, CommonData &common_data, int tag, bool jacobian=true)
OpEshelbyDynamicMaterialMomentumLhs_dX(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data, Range *forcesonlyonentities_ptr)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
OpEshelbyDynamicMaterialMomentumLhs_dv(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data, Range *forcesonlyonentities_ptr)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
OpEshelbyDynamicMaterialMomentumLhs_dx(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data, Range *forcesonlyonentities_ptr)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpEshelbyDynamicMaterialMomentumRhs(const std::string field_name, BlockData &data, CommonData &common_data, Range *forcesonlyonentities_ptr)
OpGetCommonDataAtGaussPts(const std::string field_name, CommonData &common_data)
OpGetDataAtGaussPts(const std::string field_name, std::vector< VectorDouble > &values_at_gauss_pts, std::vector< MatrixDouble > &gardient_at_gauss_pts)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
operator calculating deformation gradient
OpMassJacobian(const std::string field_name, BlockData &data, CommonData &common_data, boost::ptr_vector< MethodForForceScaling > &methods_op, int tag, bool linear=false)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpMassLhs_dM_dX(const std::string field_name, const std::string col_field, BlockData &data, CommonData &common_data)
MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
OpMassLhs_dM_dv(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data, Range *forcesonlyonentities_ptr=NULL)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
OpMassLhs_dM_dx(const std::string field_name, const std::string col_field, BlockData &data, CommonData &common_data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpMassRhs(const std::string field_name, BlockData &data, CommonData &common_data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
OpVelocityJacobian(const std::string field_name, BlockData &data, CommonData &common_data, int tag, bool jacobian=true)
OpVelocityLhs_dV_dX(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
OpVelocityLhs_dV_dv(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
OpVelocityLhs_dV_dx(const std::string vel_field, const std::string field_name, BlockData &data, CommonData &common_data)
virtual MoFEMErrorCode getJac(EntitiesFieldData::EntData &col_data, int gg)
OpVelocityRhs(const std::string field_name, BlockData &data, CommonData &common_data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
MoFEMErrorCode preProcess()
Calculate inconsistency between approximation of velocities and velocities calculated from displaceme...
UpdateAndControl(MoFEM::Interface &m_field, TS _ts, const std::string velocity_field, const std::string spatial_position_field)
MoFEMErrorCode preProcess()
Scatter values from t_u_dt on the fields.
structure grouping operators and data used for calculation of mass (convective) element \ nonlinear_e...
ConvectiveMassElement(MoFEM::Interface &m_field, short int tag)
MoFEMErrorCode setVelocityOperators(string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool ale=false)
MyVolumeFE feEnergy
calculate kinetic energy
MoFEMErrorCode setShellMatrixMassOperators(string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool linear=false)
MyVolumeFE feVelRhs
calculate right hand side for tetrahedral elements
MoFEMErrorCode addVelocityElement(string element_name, string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool ale=false, BitRefLevel bit=BitRefLevel())
MoFEMErrorCode addConvectiveMassElement(string element_name, string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool ale=false, BitRefLevel bit=BitRefLevel())
MoFEMErrorCode setKinematicEshelbyOperators(string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", Range *forces_on_entities_ptr=NULL)
MoFEMErrorCode addEshelbyDynamicMaterialMomentum(string element_name, string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool ale=false, BitRefLevel bit=BitRefLevel(), Range *intersected=NULL)
MyVolumeFE feTRhs
calculate right hand side for tetrahedral elements
MyVolumeFE feMassRhs
calculate right hand side for tetrahedral elements
MyVolumeFE feTLhs
calculate left hand side for tetrahedral elements
MoFEMErrorCode setConvectiveMassOperators(string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool ale=false, bool linear=false)
boost::ptr_vector< MethodForForceScaling > methodsOp
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
MyVolumeFE feVelLhs
calculate left hand side for tetrahedral elements
static MoFEMErrorCode applyScale(const FEMethod *fe, boost::ptr_vector< MethodForForceScaling > &methods_op, VectorDouble &nf)
Managing BitRefLevels.
Body force data structure.
virtual moab::Interface & get_moab()=0
virtual MoFEMErrorCode problem_basic_method_preProcess(const Problem *problem_ptr, BasicMethod &method, int verb=DEFAULT_VERBOSITY)=0
Set data for BasicMethod.
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
get dofs values
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1FieldData()
Return FTensor of rank 1, i.e. vector from filed data coeffinects.
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
const VectorInt & getIndices() const
Get global indices of dofs on entity.
structure to get information form mofem into EntitiesFieldData
Elastic material data structure.
intrusive_ptr for managing petsc objects
data structure for TS (time stepping) context
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Vector manager is used to create vectors \mofem_vectors.
Definition: VecManager.hpp:23