v0.14.0
Loading...
Searching...
No Matches
SpringElement.cpp
Go to the documentation of this file.
1/* \file SpringElement.cpp
2 \brief Implementation of spring boundary condition on triangular surfaces
3*/
4
5
6
7#include <MoFEM.hpp>
8using namespace MoFEM;
9
10#include <SpringElement.hpp>
11using namespace boost::numeric;
12
14 EntData &data) {
15
17 // check that the faces have associated degrees of freedom
18 const int nb_dofs = data.getIndices().size();
19 if (nb_dofs == 0)
21 if (dAta.tRis.find(getFEEntityHandle()) == dAta.tRis.end()) {
23 }
24 CHKERR commonDataPtr->getBlockData(dAta);
25
26 // size of force vector associated to the entity
27 // set equal to the number of degrees of freedom of associated with the
28 // entity
29 nF.resize(nb_dofs, false);
30 nF.clear();
31
32 // get number of Gauss points
33 const int nb_gauss_pts = data.getN().size1();
34
35 // get integration weights
37
38 // FTensor indices
42 const double normal_stiffness = commonDataPtr->springStiffnessNormal;
43 const double tangent_stiffness = commonDataPtr->springStiffnessTangent;
44
45 // create a 3d vector to be used as the normal to the face with length equal
46 // to the face area
48 FTensor::Tensor2<double, 3, 3> t_normal_projection;
49 FTensor::Tensor2<double, 3, 3> t_tangent_projection;
50
51 // Extract solution at Gauss points
52 auto t_solution_at_gauss_point =
53 getFTensor1FromMat<3>(*commonDataPtr->xAtPts);
54 auto t_init_solution_at_gauss_point =
55 getFTensor1FromMat<3>(*commonDataPtr->xInitAtPts);
56 FTensor::Tensor1<double, 3> t_displacement_at_gauss_point;
57
58 auto t_normal_1 = getFTensor1FromMat<3>(commonDataPtr->normalVector);
59 // loop over all Gauss points of the face
60 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
61 t_normal(i) = t_normal_1(i);
62
63 const double normal_length = std::sqrt(t_normal(k) * t_normal(k));
64 t_normal(i) = t_normal(i) / normal_length;
65 t_normal_projection(i, j) = t_normal(i) * t_normal(j);
66 t_tangent_projection(i, j) = t_normal_projection(i, j);
67 t_tangent_projection(0, 0) -= 1;
68 t_tangent_projection(1, 1) -= 1;
69 t_tangent_projection(2, 2) -= 1;
70 // Calculate the displacement at the Gauss point
71 if (is_spatial_position) { // "SPATIAL_POSITION"
72 t_displacement_at_gauss_point(i) =
73 t_solution_at_gauss_point(i) - t_init_solution_at_gauss_point(i);
74 } else { // e.g. "DISPLACEMENT" or "U"
75 t_displacement_at_gauss_point(i) = t_solution_at_gauss_point(i);
76 }
77
78 double w = t_w * normal_length; // area was constant
79 if (getNumeredEntFiniteElementPtr()->getEntType() == MBTRI) {
80 w *= 0.5;
81 }
82
83 auto t_base_func = data.getFTensor0N(gg, 0);
84
85 // create a vector t_nf whose pointer points an array of 3 pointers
86 // pointing to nF memory location of components
88 &nF[2]);
89
90 for (int rr = 0; rr != nb_dofs / 3; ++rr) { // loop over the nodes
91
92 t_nf(i) += w * t_base_func * normal_stiffness *
93 t_normal_projection(i, j) * t_displacement_at_gauss_point(j);
94 t_nf(i) -= w * t_base_func * tangent_stiffness *
95 t_tangent_projection(i, j) * t_displacement_at_gauss_point(j);
96
97 // move to next base function
98 ++t_base_func;
99 // move the pointer to next element of t_nf
100 ++t_nf;
101 }
102 // move to next integration weight
103 ++t_w;
104 // move to the solutions at the next Gauss point
105 ++t_solution_at_gauss_point;
106 ++t_init_solution_at_gauss_point;
107 ++t_normal_1;
108 }
109 // add computed values of spring in the global right hand side vector
110 Vec f = getFEMethod()->ksp_f != PETSC_NULL ? getFEMethod()->ksp_f
111 : getFEMethod()->snes_f;
112 CHKERR VecSetValues(f, nb_dofs, &data.getIndices()[0], &nF[0], ADD_VALUES);
114}
115
117 EntityType row_type,
118 EntityType col_type,
119 EntData &row_data,
120 EntData &col_data) {
122 // check if the volumes have associated degrees of freedom
123 const int row_nb_dofs = row_data.getIndices().size();
124 if (!row_nb_dofs)
126
127 const int col_nb_dofs = col_data.getIndices().size();
128 if (!col_nb_dofs)
130
131 if (dAta.tRis.find(getFEEntityHandle()) == dAta.tRis.end()) {
133 }
134
135 CHKERR commonDataPtr->getBlockData(dAta);
136 // size associated to the entity
137 locKs.resize(row_nb_dofs, col_nb_dofs, false);
138 locKs.clear();
139
140 // get number of Gauss points
141 const int row_nb_gauss_pts = row_data.getN().size1();
142 if (!row_nb_gauss_pts) // check if number of Gauss point <> 0
144
145 // get intergration weights
146 auto t_w = getFTensor0IntegrationWeight();
147
151
152 // create a 3d vector to be used as the normal to the face with length equal
153 // to the face area
155
156 // First tangent vector
157 auto t_tangent1_ptr = getFTensor1Tangent1AtGaussPts();
159 t_tangent1(i) = t_tangent1_ptr(i);
160
161 // Second tangent vector, such that t_n = t_t1 x t_t2 | t_t2 = t_n x t_t1
163 t_tangent2(i) = FTensor::levi_civita(i, j, k) * t_normal(j) * t_tangent1(k);
164
165 auto t_tangent2_ptr = getFTensor1Tangent2AtGaussPts();
166 t_normal(i) =
167 FTensor::levi_civita(i, j, k) * t_tangent1(j) * t_tangent2_ptr(k);
168
169 auto normal_original_slave =
170 getFTensor1FromMat<3>(commonDataPtr->normalVector);
171
172 FTensor::Tensor2<double, 3, 3> t_normal_projection;
173 FTensor::Tensor2<double, 3, 3> t_tangent_projection;
174 // loop over the Gauss points
175 for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
176 t_normal(i) = normal_original_slave(i);
177 const double normal_length = std::sqrt(t_normal(k) * t_normal(k));
178 t_normal(i) = t_normal(i) / normal_length;
179 t_normal_projection(i, j) = t_normal(i) * t_normal(j);
180 t_tangent_projection(i, j) = t_normal_projection(i, j);
181 t_tangent_projection(0, 0) -= 1;
182 t_tangent_projection(1, 1) -= 1;
183 t_tangent_projection(2, 2) -= 1;
184
185 // get area and integration weight
186 double w = t_w * normal_length;
187 if (getNumeredEntFiniteElementPtr()->getEntType() == MBTRI) {
188 w *= 0.5;
189 }
190
191 auto t_row_base_func = row_data.getFTensor0N(gg, 0);
192
193 for (int rr = 0; rr != row_nb_dofs / 3; rr++) {
194 auto assemble_m = getFTensor2FromArray<3, 3, 3>(locKs, 3 * rr);
195 auto t_col_base_func = col_data.getFTensor0N(gg, 0);
196 for (int cc = 0; cc != col_nb_dofs / 3; cc++) {
197 assemble_m(i, j) += w * t_row_base_func * t_col_base_func *
198 commonDataPtr->springStiffnessNormal *
199 t_normal_projection(i, j);
200 assemble_m(i, j) -= w * t_row_base_func * t_col_base_func *
201 commonDataPtr->springStiffnessTangent *
202 t_tangent_projection(i, j);
203 ++t_col_base_func;
204 ++assemble_m;
205 }
206 ++t_row_base_func;
207 }
208 // move to next integration weight
209 ++t_w;
210 ++normal_original_slave;
211 }
212
213 // Add computed values of spring stiffness to the global LHS matrix
214 CHKERR MatSetValues(getKSPB(), row_data, col_data, &locKs(0, 0), ADD_VALUES);
215
216 // is symmetric
217 if (row_side != col_side || row_type != col_type) {
218 transLocKs.resize(col_nb_dofs, row_nb_dofs, false);
219 noalias(transLocKs) = trans(locKs);
220
221 CHKERR MatSetValues(getKSPB(), col_data, row_data, &transLocKs(0, 0),
222 ADD_VALUES);
223 }
224
226}
227
229 EntityType row_type,
230 EntityType col_type,
231 EntData &row_data,
232 EntData &col_data) {
234
235 // check if the volumes have associated degrees of freedom
236 const int row_nb_dofs = row_data.getIndices().size();
237 if (!row_nb_dofs)
239
240 const int col_nb_dofs = col_data.getIndices().size();
241 if (!col_nb_dofs)
243
244 if (dAta.tRis.find(getFEEntityHandle()) == dAta.tRis.end()) {
246 }
247
248 CHKERR commonDataPtr->getBlockData(dAta);
249 // size associated to the entity
250 locKs.resize(row_nb_dofs, col_nb_dofs, false);
251 locKs.clear();
252
253 // get number of Gauss points
254 const int row_nb_gauss_pts = row_data.getN().size1();
255 if (!row_nb_gauss_pts) // check if number of Gauss point <> 0
257
258 // get intergration weights
259 auto t_w = getFTensor0IntegrationWeight();
260
265
266 auto make_vec_der = [&](auto t_N, auto t_1, auto t_2) {
268 t_n(i, j) = 0;
269 t_n(i, j) += FTensor::levi_civita(i, j, k) * t_2(k) * t_N(0);
270 t_n(i, j) -= FTensor::levi_civita(i, j, k) * t_1(k) * t_N(1);
271 return t_n;
272 };
273
274 auto make_vec_der_2 = [&](auto t_N, auto t_1, auto t_2) {
276 t_normal(i) = FTensor::levi_civita(i, j, k) * t_1(j) * t_2(k);
277 const double normal_norm = std::sqrt(t_normal(i) * t_normal(i));
279 t_n(i, j) = 0;
280 t_n(i, j) += FTensor::levi_civita(i, j, k) * t_2(k) * t_N(0);
281 t_n(i, j) -= FTensor::levi_civita(i, j, k) * t_1(k) * t_N(1);
283 t_result(i, j) = 0;
284 t_result(i, j) =
285 t_n(i, j) / normal_norm - t_normal(i) * t_n(k, j) * t_normal(k) /
286 (normal_norm * normal_norm * normal_norm);
287 return t_result;
288 };
289
290 const double normal_stiffness = commonDataPtr->springStiffnessNormal;
291 const double tangent_stiffness = commonDataPtr->springStiffnessTangent;
292
293 // create a 3d vector to be used as the normal to the face with length equal
294 // to the face area
296 FTensor::Tensor2<double, 3, 3> t_normal_projection;
297 FTensor::Tensor2<double, 3, 3> t_tangent_projection;
298 FTensor::Tensor1<double, 3> t_displacement_at_gauss_point;
299
300 // Extract solution at Gauss points
301 auto t_solution_at_gauss_point =
302 getFTensor1FromMat<3>(*commonDataPtr->xAtPts);
303 auto t_init_solution_at_gauss_point =
304 getFTensor1FromMat<3>(*commonDataPtr->xInitAtPts);
305
306 auto t_1 = getFTensor1FromMat<3>(commonDataPtr->tangent1);
307 auto t_2 = getFTensor1FromMat<3>(commonDataPtr->tangent2);
308
309 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
310
311 auto normal_at_gp = getFTensor1FromMat<3>(commonDataPtr->normalVector);
312
313 // loop over the Gauss points
314 for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
315 // get area and integration weight
316 t_normal(i) = normal_at_gp(i);
317 const double normal_length = std::sqrt(t_normal(k) * t_normal(k));
318 t_normal(i) = t_normal(i) / normal_length;
319
320 t_normal_projection(i, j) = t_normal(i) * t_normal(j);
321 t_tangent_projection(i, j) = t_normal_projection(i, j);
322 t_tangent_projection(0, 0) -= 1;
323 t_tangent_projection(1, 1) -= 1;
324 t_tangent_projection(2, 2) -= 1;
325
326 double w = 0.5 * t_w;
327
328 t_displacement_at_gauss_point(i) =
329 t_solution_at_gauss_point(i) - t_init_solution_at_gauss_point(i);
330
331 auto t_row_base_func = row_data.getFTensor0N(gg, 0);
332
333 for (int rr = 0; rr != row_nb_dofs / 3; rr++) {
334 auto t_col_base_func = col_data.getFTensor0N(gg, 0);
335 auto t_N = col_data.getFTensor1DiffN<2>(gg, 0);
336 auto assemble_m = getFTensor2FromArray<3, 3, 3>(locKs, 3 * rr);
337 for (int cc = 0; cc != col_nb_dofs / 3; cc++) {
338 auto t_d_n = make_vec_der(t_N, t_1, t_2);
339 auto t_d_n_2 = make_vec_der_2(t_N, t_1, t_2);
340
341 assemble_m(i, j) -= w * normal_length * t_col_base_func *
342 normal_stiffness * t_row_base_func *
343 t_normal_projection(i, j);
344
345 assemble_m(i, j) += w * normal_length * t_col_base_func *
346 tangent_stiffness * t_row_base_func *
347 t_tangent_projection(i, j);
348
349 assemble_m(i, j) +=
350 w * (normal_stiffness - tangent_stiffness) * t_row_base_func *
351 (t_d_n(i, j) * (t_normal(k) * t_displacement_at_gauss_point(k)) +
352 normal_length * t_normal(i) *
353 (t_d_n_2(k, j) * t_displacement_at_gauss_point(k)));
354 // constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
355 assemble_m(i, j) += w * (t_d_n(l, j) * t_normal(l)) *
356 tangent_stiffness * t_row_base_func * t_kd(i, k) *
357 t_displacement_at_gauss_point(k);
358
359 ++t_col_base_func;
360 ++t_N;
361 ++assemble_m;
362 }
363 ++t_row_base_func;
364 }
365 // move to next integration weight
366 ++t_w;
367 ++t_solution_at_gauss_point;
368 ++t_init_solution_at_gauss_point;
369 ++t_1;
370 ++t_2;
371 ++normal_at_gp;
372 }
373
374 // Add computed values of spring stiffness to the global LHS matrix
375 CHKERR MatSetValues(getKSPB(), row_data, col_data, &locKs(0, 0), ADD_VALUES);
376
378}
379
381 int row_side, int col_side, EntityType row_type, EntityType col_type,
382 EntData &row_data, EntData &col_data) {
383
385
386 if (dataAtSpringPts->faceRowData == nullptr)
388
389 if (row_type != MBVERTEX)
391
392 row_nb_dofs = dataAtSpringPts->faceRowData->getIndices().size();
393 if (!row_nb_dofs)
395 col_nb_dofs = col_data.getIndices().size();
396 if (!col_nb_dofs)
398
399 nb_gauss_pts = dataAtSpringPts->faceRowData->getN().size1();
400
401 nb_base_fun_row = dataAtSpringPts->faceRowData->getFieldData().size() / 3;
402 nb_base_fun_col = col_data.getFieldData().size() / 3;
403
404 NN.resize(3 * nb_base_fun_row, 3 * nb_base_fun_col, false);
405 NN.clear();
406
407 // integrate local matrix for entity block
408 CHKERR iNtegrate(*(dataAtSpringPts->faceRowData), col_data);
409
410 // assemble local matrix
411 CHKERR aSsemble(*(dataAtSpringPts->faceRowData), col_data);
412
414}
415
418 EntData &col_data) {
419
421
422 // get pointer to first global index on row
423 const int *row_indices = &*row_data.getIndices().data().begin();
424 // get pointer to first global index on column
425 const int *col_indices = &*col_data.getIndices().data().begin();
426
427 auto &data = *dataAtSpringPts;
428
429 rowIndices.resize(row_nb_dofs, false);
430 noalias(rowIndices) = row_data.getIndices();
431 row_indices = &rowIndices[0];
432 VectorDofs &dofs = row_data.getFieldDofs();
433 VectorDofs::iterator dit = dofs.begin();
434 for (int ii = 0; dit != dofs.end(); ++dit, ++ii) {
435 if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
436 data.forcesOnlyOnEntitiesRow.end()) {
437 rowIndices[ii] = -1;
438 }
439 }
440
441 Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
442 : getFEMethod()->snes_B;
443 // assemble local matrix
444 CHKERR MatSetValues(B, row_nb_dofs, row_indices, col_nb_dofs, col_indices,
445 &*NN.data().begin(), ADD_VALUES);
446
448}
449
451 EntData &row_data, EntData &col_data) {
452
454
459
460 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
462 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
463 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
464 &m(r + 2, c + 2));
465 };
467 FTensor::Tensor2<double, 3, 3> t_normal_projection;
468 FTensor::Tensor2<double, 3, 3> t_tangent_projection;
469
470 const double normal_stiffness = dataAtSpringPts->springStiffnessNormal;
471 const double tangent_stiffness = dataAtSpringPts->springStiffnessTangent;
472
473 // Extract solution at Gauss points
474 auto t_solution_at_gauss_point =
475 getFTensor1FromMat<3>(*dataAtSpringPts->xAtPts);
476 auto t_init_solution_at_gauss_point =
477 getFTensor1FromMat<3>(*dataAtSpringPts->xInitAtPts);
478 FTensor::Tensor1<double, 3> t_displacement_at_gauss_point;
479
480 auto t_normal_1 = getFTensor1FromMat<3>(dataAtSpringPts->normalVector);
481 auto t_w = getFTensor0IntegrationWeight();
482
483 auto t_inv_H = getFTensor2FromMat<3, 3>(*dataAtSpringPts->invHMat);
484 auto t_F = getFTensor2FromMat<3, 3>(*dataAtSpringPts->FMat);
485 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
486
487 t_normal(i) = t_normal_1(i);
488
489 const double normal_length = std::sqrt(t_normal(k) * t_normal(k));
490 t_normal(i) = t_normal(i) / normal_length;
491 t_normal_projection(i, j) = t_normal(i) * t_normal(j);
492 t_tangent_projection(i, j) = t_normal_projection(i, j);
493 t_tangent_projection(0, 0) -= 1;
494 t_tangent_projection(1, 1) -= 1;
495 t_tangent_projection(2, 2) -= 1;
496
497 t_displacement_at_gauss_point(i) =
498 t_solution_at_gauss_point(i) - t_init_solution_at_gauss_point(i);
499
500 double w = 0.5 * t_w * normal_length;
501
502 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
503 FTensor::Tensor0<double *> t_col_base(&col_data.getN()(gg, 0));
504 int bbc = 0;
505 for (; bbc != nb_base_fun_col; ++bbc) {
506
507 FTensor::Tensor0<double *> t_row_base(&row_data.getN()(gg, 0));
508
509 int bbr = 0;
510 for (; bbr != nb_base_fun_row; ++bbr) {
511
512 auto t_assemble = get_tensor2(NN, 3 * bbr, 3 * bbc);
513
514 t_assemble(i, j) -= w * t_row_base * t_col_base * normal_stiffness *
515 t_F(k, i) * t_normal_projection(k, j);
516 t_assemble(i, j) += w * t_row_base * t_col_base * tangent_stiffness *
517 t_F(k, i) * t_tangent_projection(k, j);
518
519 t_assemble(i, j) -= w * normal_stiffness * t_row_base * t_inv_H(k, i) *
520 t_col_diff_base(k) * t_normal_projection(j, l) *
521 t_displacement_at_gauss_point(l);
522 t_assemble(i, j) += w * tangent_stiffness * t_row_base * t_inv_H(k, i) *
523 t_col_diff_base(k) * t_tangent_projection(j, l) *
524 t_displacement_at_gauss_point(l);
525
526 ++t_row_base;
527 }
528 ++t_col_diff_base;
529 ++t_col_base;
530 }
531 ++t_w;
532 ++t_solution_at_gauss_point;
533 ++t_init_solution_at_gauss_point;
534 ++t_normal_1;
535 ++t_inv_H;
536 ++t_F;
537 }
538
540}
541
544 EntData &col_data) {
545
547
548 // get pointer to first global index on row
549 const int *row_indices = &*row_data.getIndices().data().begin();
550 // get pointer to first global index on column
551 const int *col_indices = &*col_data.getIndices().data().begin();
552
553 auto &data = *dataAtSpringPts;
554
555 rowIndices.resize(row_nb_dofs, false);
556 noalias(rowIndices) = row_data.getIndices();
557 row_indices = &rowIndices[0];
558 VectorDofs &dofs = row_data.getFieldDofs();
559 VectorDofs::iterator dit = dofs.begin();
560 for (int ii = 0; dit != dofs.end(); ++dit, ++ii) {
561 if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
562 data.forcesOnlyOnEntitiesRow.end()) {
563 rowIndices[ii] = -1;
564 }
565 }
566
567 Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
568 : getFEMethod()->snes_B;
569 // assemble local matrix
570 CHKERR MatSetValues(B, row_nb_dofs, row_indices, col_nb_dofs, col_indices,
571 &*NN.data().begin(), ADD_VALUES);
572
574}
575
577 int row_side, int col_side, EntityType row_type, EntityType col_type,
578 EntData &row_data, EntData &col_data) {
579
581
582 row_nb_dofs = row_data.getIndices().size();
583 if (!row_nb_dofs)
585 col_nb_dofs = col_data.getIndices().size();
586 if (!col_nb_dofs)
588 nb_gauss_pts = row_data.getN().size1();
589
590 nb_base_fun_row = row_data.getFieldData().size() / 3;
591 nb_base_fun_col = col_data.getFieldData().size() / 3;
592
593 NN.resize(3 * nb_base_fun_row, 3 * nb_base_fun_col, false);
594 NN.clear();
595
596 diagonal_block = (row_type == col_type) && (row_side == col_side);
597
598 if (col_type == MBVERTEX) {
599 dataAtSpringPts->faceRowData = &row_data;
600 CHKERR loopSideVolumes(sideFeName, *sideFe);
601 }
602
603 // integrate local matrix for entity block
604 CHKERR iNtegrate(row_data, col_data);
605
606 // assemble local matrix
607 CHKERR aSsemble(row_data, col_data);
608
610}
611
614 EntData &col_data) {
615
617
622
623 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
625 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
626 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
627 &m(r + 2, c + 2));
628 };
629
630 auto make_vec_der = [&](auto t_N, auto t_1, auto t_2) {
632 t_n(i, j) = 0;
633 t_n(i, j) += FTensor::levi_civita(i, j, k) * t_2(k) * t_N(0);
634 t_n(i, j) -= FTensor::levi_civita(i, j, k) * t_1(k) * t_N(1);
635 return t_n;
636 };
637
638 auto make_vec_der_2 = [&](auto t_N, auto t_1, auto t_2) {
640 t_normal(i) = FTensor::levi_civita(i, j, k) * t_1(j) * t_2(k);
641 const double normal_norm = std::sqrt(t_normal(i) * t_normal(i));
643 t_n(i, j) = 0;
644 t_n(i, j) += FTensor::levi_civita(i, j, k) * t_2(k) * t_N(0);
645 t_n(i, j) -= FTensor::levi_civita(i, j, k) * t_1(k) * t_N(1);
647 t_result(i, j) = 0;
648 t_result(i, j) =
649 t_n(i, j) / normal_norm - t_normal(i) * t_n(k, j) * t_normal(k) /
650 (normal_norm * normal_norm * normal_norm);
651 return t_result;
652 };
653
654 dataAtSpringPts->faceRowData = nullptr;
655 CHKERR loopSideVolumes(sideFeName, *sideFe);
656
657 auto t_F = getFTensor2FromMat<3, 3>(*dataAtSpringPts->FMat);
658
659 auto t_w = getFTensor0IntegrationWeight();
660 auto t_1 = getFTensor1FromMat<3>(dataAtSpringPts->tangent1);
661 auto t_2 = getFTensor1FromMat<3>(dataAtSpringPts->tangent2);
662
663 const double normal_stiffness = dataAtSpringPts->springStiffnessNormal;
664 const double tangent_stiffness = dataAtSpringPts->springStiffnessTangent;
665
666 // create a 3d vector to be used as the normal to the face with length equal
667 // to the face area
669 FTensor::Tensor2<double, 3, 3> t_normal_projection;
670 FTensor::Tensor2<double, 3, 3> t_tangent_projection;
671
672 // Extract solution at Gauss points
673 auto t_solution_at_gauss_point =
674 getFTensor1FromMat<3>(*dataAtSpringPts->xAtPts);
675 auto t_init_solution_at_gauss_point =
676 getFTensor1FromMat<3>(*dataAtSpringPts->xInitAtPts);
677 FTensor::Tensor1<double, 3> t_displacement_at_gauss_point;
678
679 auto t_normal_1 = getFTensor1FromMat<3>(dataAtSpringPts->normalVector);
681
682 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
683 t_normal(i) = t_normal_1(i);
684
685 const double normal_length = std::sqrt(t_normal(k) * t_normal(k));
686 t_normal(i) = t_normal(i) / normal_length;
687 t_normal_projection(i, j) = t_normal(i) * t_normal(j);
688 t_tangent_projection(i, j) = t_normal_projection(i, j);
689 t_tangent_projection(0, 0) -= 1;
690 t_tangent_projection(1, 1) -= 1;
691 t_tangent_projection(2, 2) -= 1;
692 // Calculate the displacement at the Gauss point
693 t_displacement_at_gauss_point(i) =
694 t_solution_at_gauss_point(i) - t_init_solution_at_gauss_point(i);
695 double val = 0.5 * t_w;
696
697 auto t_N = col_data.getFTensor1DiffN<2>(gg, 0);
698 auto t_col_base_func = col_data.getFTensor0N(gg, 0);
699 int bbc = 0;
700 for (; bbc != nb_base_fun_col; ++bbc) {
701
702 auto t_row_base_func = row_data.getFTensor0N(gg, 0);
703
704 int bbr = 0;
705 for (; bbr != nb_base_fun_row; ++bbr) {
706
707 auto t_d_n = make_vec_der(t_N, t_1, t_2);
708
709 auto t_assemble = get_tensor2(NN, 3 * bbr, 3 * bbc);
710
711 auto t_d_n_2 = make_vec_der_2(t_N, t_1, t_2);
712
713 t_assemble(i, j) += val * normal_length * t_col_base_func *
714 normal_stiffness * t_row_base_func * t_F(k, i) *
715 t_normal_projection(k, j);
716
717 t_assemble(i, j) -= val * normal_length * t_col_base_func *
718 tangent_stiffness * t_row_base_func * t_F(k, i) *
719 t_tangent_projection(k, j);
720
721 t_assemble(i, j) -=
722 val * (normal_stiffness - tangent_stiffness) * t_row_base_func *
723 t_F(l, i) *
724 (t_d_n(l, j) * (t_normal(k) * t_displacement_at_gauss_point(k)) +
725 normal_length * t_normal(l) *
726 (t_d_n_2(k, j) * t_displacement_at_gauss_point(k)));
727 // constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
728 t_assemble(i, j) -= val * (t_d_n(l, j) * t_normal(l)) *
729 tangent_stiffness * t_row_base_func * t_F(k, i) *
730 t_displacement_at_gauss_point(k);
731
732 ++t_row_base_func;
733 }
734 ++t_N;
735 ++t_col_base_func;
736 }
737 ++t_F;
738 ++t_w;
739 ++t_1;
740 ++t_2;
741 ++t_normal_1;
742 ++t_solution_at_gauss_point;
743 ++t_init_solution_at_gauss_point;
744 }
745
747}
748
750 EntData &row_data, EntData &col_data) {
751
753
760
761 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
763 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
764 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
765 &m(r + 2, c + 2));
766 };
767
768 const double normal_stiffness = dataAtSpringPts->springStiffnessNormal;
769 const double tangent_stiffness = dataAtSpringPts->springStiffnessTangent;
770
771 // create a 3d vector to be used as the normal to the face with length equal
772 // to the face area
774 FTensor::Tensor2<double, 3, 3> t_normal_projection;
775 FTensor::Tensor2<double, 3, 3> t_tangent_projection;
776
777 // Extract solution at Gauss points
778 auto t_solution_at_gauss_point =
779 getFTensor1FromMat<3>(*dataAtSpringPts->xAtPts);
780 auto t_init_solution_at_gauss_point =
781 getFTensor1FromMat<3>(*dataAtSpringPts->xInitAtPts);
782 FTensor::Tensor1<double, 3> t_displacement_at_gauss_point;
783
784 auto t_normal_1 = getFTensor1FromMat<3>(dataAtSpringPts->normalVector);
785
786 auto t_w = getFTensor0IntegrationWeight();
787
788 auto t_h = getFTensor2FromMat<3, 3>(*dataAtSpringPts->hMat);
789 auto t_inv_H = getFTensor2FromMat<3, 3>(*dataAtSpringPts->invHMat);
790
792
793 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
794 t_normal(i) = t_normal_1(i);
795
796 const double normal_length = std::sqrt(t_normal(k) * t_normal(k));
797 t_normal(i) = t_normal(i) / normal_length;
798 t_normal_projection(i, j) = t_normal(i) * t_normal(j);
799 t_tangent_projection(i, j) = t_normal_projection(i, j);
800 t_tangent_projection(0, 0) -= 1;
801 t_tangent_projection(1, 1) -= 1;
802 t_tangent_projection(2, 2) -= 1;
803 // Calculate the displacement at the Gauss point
804 t_displacement_at_gauss_point(i) =
805 t_solution_at_gauss_point(i) - t_init_solution_at_gauss_point(i);
806
807 double val = 0.5 * t_w * normal_length;
808
809 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
810
811 int bbc = 0;
812 for (; bbc != nb_base_fun_col; ++bbc) {
813
814 FTensor::Tensor0<double *> t_row_base(&row_data.getN()(gg, 0));
815
816 int bbr = 0;
817 for (; bbr != nb_base_fun_row; ++bbr) {
818
819 auto t_assemble = get_tensor2(NN, 3 * bbr, 3 * bbc);
820
821 t_assemble(i, j) +=
822 val * t_row_base * normal_stiffness * t_inv_H(l, j) *
823 t_col_diff_base(m) * t_inv_H(m, i) * t_h(k, l) *
824 (t_normal_projection(k, q) * t_displacement_at_gauss_point(q));
825
826 t_assemble(i, j) -=
827 val * t_row_base * tangent_stiffness * t_inv_H(l, j) *
828 t_col_diff_base(m) * t_inv_H(m, i) * t_h(k, l) *
829 (t_tangent_projection(k, q) * t_displacement_at_gauss_point(q));
830
831 ++t_row_base;
832 }
833 ++t_col_diff_base;
834 }
835 ++t_w;
836 ++t_h;
837 ++t_inv_H;
838 ++t_solution_at_gauss_point;
839 ++t_init_solution_at_gauss_point;
840 ++t_normal_1;
841 }
842
844}
845
847 EntityType type,
848 EntData &row_data) {
849
854
855 // get number of integration points
856 const int nb_integration_pts = getGaussPts().size2();
857
858 auto t_h = getFTensor2FromMat<3, 3>(*commonDataPtr->hMat);
859 auto t_H = getFTensor2FromMat<3, 3>(*commonDataPtr->HMat);
860
861 commonDataPtr->detHVec->resize(nb_integration_pts, false);
862 commonDataPtr->invHMat->resize(9, nb_integration_pts, false);
863 commonDataPtr->FMat->resize(9, nb_integration_pts, false);
864
865 commonDataPtr->detHVec->clear();
866 commonDataPtr->invHMat->clear();
867 commonDataPtr->FMat->clear();
868
869 auto t_detH = getFTensor0FromVec(*commonDataPtr->detHVec);
870 auto t_invH = getFTensor2FromMat<3, 3>(*commonDataPtr->invHMat);
871 auto t_F = getFTensor2FromMat<3, 3>(*commonDataPtr->FMat);
872
873 for (int gg = 0; gg != nb_integration_pts; ++gg) {
874 CHKERR determinantTensor3by3(t_H, t_detH);
875 CHKERR invertTensor3by3(t_H, t_detH, t_invH);
876 t_F(i, j) = t_h(i, k) * t_invH(k, j);
877 ++t_h;
878 ++t_H;
879 ++t_detH;
880 ++t_invH;
881 ++t_F;
882 }
883
885}
886
888 EntityType type,
889 EntData &row_data) {
890
892
893 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
894 dAta.tRis.end()) {
896 }
897
898 // get number of dofs on row
899 nbRows = row_data.getIndices().size();
900 // if no dofs on row, exit that work, nothing to do here
901 if (!nbRows)
903
904 nF.resize(nbRows, false);
905 nF.clear();
906
907 // get number of integration points
908 nbIntegrationPts = getGaussPts().size2();
909
910 // integrate local matrix for entity block
911 CHKERR iNtegrate(row_data);
912
913 // assemble local matrix
914 CHKERR aSsemble(row_data);
915
917}
918
921
922 CHKERR loopSideVolumes(sideFeName, *sideFe);
923
924 // check that the faces have associated degrees of freedom
925 const int nb_dofs = data.getIndices().size();
926 if (nb_dofs == 0)
928
929 CHKERR dataAtPts->getBlockData(dAta);
930
931 // get integration weights
932 auto t_w = getFTensor0IntegrationWeight();
933
934 // FTensor indices
938
939 const double normal_stiffness = dataAtPts->springStiffnessNormal;
940 const double tangent_stiffness = dataAtPts->springStiffnessTangent;
941
942 // create a 3d vector to be used as the normal to the face with length equal
943 // to the face area
945 FTensor::Tensor2<double, 3, 3> t_normal_projection;
946 FTensor::Tensor2<double, 3, 3> t_tangent_projection;
947
948 // Extract solution at Gauss points
949 auto t_solution_at_gauss_point = getFTensor1FromMat<3>(*dataAtPts->xAtPts);
950 auto t_init_solution_at_gauss_point =
951 getFTensor1FromMat<3>(*dataAtPts->xInitAtPts);
952 FTensor::Tensor1<double, 3> t_displacement_at_gauss_point;
953
954 auto t_normal_1 = getFTensor1FromMat<3>(dataAtPts->normalVector);
955
956 auto t_F = getFTensor2FromMat<3, 3>(*dataAtPts->FMat);
957
958 // loop over all Gauss points of the face
959 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
960 t_normal(i) = t_normal_1(i);
961
962 const double normal_length = std::sqrt(t_normal(k) * t_normal(k));
963 t_normal(i) = t_normal(i) / normal_length;
964 t_normal_projection(i, j) = t_normal(i) * t_normal(j);
965 t_tangent_projection(i, j) = t_normal_projection(i, j);
966 t_tangent_projection(0, 0) -= 1;
967 t_tangent_projection(1, 1) -= 1;
968 t_tangent_projection(2, 2) -= 1;
969 // Calculate the displacement at the Gauss point
970 t_displacement_at_gauss_point(i) =
971 t_solution_at_gauss_point(i) - t_init_solution_at_gauss_point(i);
972
973 double w = t_w * 0.5 * normal_length; // area was constant
974
975 auto t_base_func = data.getFTensor0N(gg, 0);
976
977 // create a vector t_nf whose pointer points an array of 3 pointers
978 // pointing to nF memory location of components
980 &nF[2]);
981
982 for (int rr = 0; rr != nb_dofs / 3; ++rr) { // loop over the nodes
983
984 t_nf(i) -= w * t_base_func * normal_stiffness * t_F(k, i) *
985 t_normal_projection(k, j) * t_displacement_at_gauss_point(j);
986 t_nf(i) += w * t_base_func * tangent_stiffness * t_F(k, i) *
987 t_tangent_projection(k, j) * t_displacement_at_gauss_point(j);
988
989 // move to next base function
990 ++t_base_func;
991 // move the pointer to next element of t_nf
992 ++t_nf;
993 }
994 // move to next integration weight
995 ++t_w;
996 // move to the solutions at the next Gauss point
997 ++t_solution_at_gauss_point;
998 ++t_init_solution_at_gauss_point;
999 ++t_normal_1;
1000 ++t_F;
1001 }
1002 // add computed values of spring in the global right hand side vector
1003 // Vec f = getFEMethod()->ksp_f != PETSC_NULL ? getFEMethod()->ksp_f
1004 // : getFEMethod()->snes_f;
1005 // CHKERR VecSetValues(f, nb_dofs, &data.getIndices()[0], &nF[0], ADD_VALUES);
1007}
1008
1011
1012 // get pointer to first global index on row
1013 const int *row_indices = &*row_data.getIndices().data().begin();
1014
1015 auto &data = *dataAtPts;
1016
1017 rowIndices.resize(nbRows, false);
1018 noalias(rowIndices) = row_data.getIndices();
1019 row_indices = &rowIndices[0];
1020 VectorDofs &dofs = row_data.getFieldDofs();
1021 VectorDofs::iterator dit = dofs.begin();
1022 for (int ii = 0; dit != dofs.end(); ++dit, ++ii) {
1023 if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
1024 data.forcesOnlyOnEntitiesRow.end()) {
1025 rowIndices[ii] = -1;
1026 }
1027 }
1028
1029 CHKERR VecSetOption(getSNESf(), VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
1030 CHKERR VecSetValues(getSNESf(), nbRows, row_indices, &*nF.data().begin(),
1031 ADD_VALUES);
1032
1034}
1035
1037 EntityType type,
1038 EntData &data) {
1040
1041 if (data.getFieldData().size() == 0)
1042 PetscFunctionReturn(0);
1043
1044 ngp = data.getN().size1();
1045
1046 unsigned int nb_dofs = data.getFieldData().size() / 3;
1048 if (type == MBVERTEX) {
1049 dataAtIntegrationPts->tangent1.resize(3, ngp, false);
1050 dataAtIntegrationPts->tangent1.clear();
1051
1052 dataAtIntegrationPts->tangent2.resize(3, ngp, false);
1053 dataAtIntegrationPts->tangent2.clear();
1054 }
1055
1056 auto t_1 = getFTensor1FromMat<3>(dataAtIntegrationPts->tangent1);
1057 auto t_2 = getFTensor1FromMat<3>(dataAtIntegrationPts->tangent2);
1058
1059 for (unsigned int gg = 0; gg != ngp; ++gg) {
1060 auto t_N = data.getFTensor1DiffN<2>(gg, 0);
1061 auto t_dof = data.getFTensor1FieldData<3>();
1062
1063 for (unsigned int dd = 0; dd != nb_dofs; ++dd) {
1064 t_1(i) += t_dof(i) * t_N(0);
1065 t_2(i) += t_dof(i) * t_N(1);
1066 ++t_dof;
1067 ++t_N;
1068 }
1069 ++t_1;
1070 ++t_2;
1071 }
1072
1074}
1075
1077 EntData &data) {
1079
1080 if (data.getFieldData().size() == 0)
1082
1083 ngp = data.getN().size1();
1084
1088
1089 if (type == MBVERTEX) {
1090 dataAtIntegrationPts->normalVector.resize(3, ngp, false);
1091 dataAtIntegrationPts->normalVector.clear();
1092 }
1093
1094 auto normal_original_slave =
1095 getFTensor1FromMat<3>(dataAtIntegrationPts->normalVector);
1096
1097 auto tangent_1 = getFTensor1FromMat<3>(dataAtIntegrationPts->tangent1);
1098 auto tangent_2 = getFTensor1FromMat<3>(dataAtIntegrationPts->tangent2);
1099
1100 for (unsigned int gg = 0; gg != ngp; ++gg) {
1101 normal_original_slave(i) =
1102 FTensor::levi_civita(i, j, k) * tangent_1(j) * tangent_2(k);
1103 ++normal_original_slave;
1104 ++tangent_1;
1105 ++tangent_2;
1106 }
1107
1109}
1110
1112 int row_side, int col_side, EntityType row_type, EntityType col_type,
1113 EntData &row_data, EntData &col_data) {
1114
1116
1117 if (col_type != MBVERTEX)
1119
1120 dataAtSpringPts->faceRowData = &row_data;
1121 CHKERR loopSideVolumes(sideFeName, *sideFe);
1122
1124}
1125
1128 const std::string field_name,
1129 const std::string mesh_nodals_positions) {
1131
1132 // Define boundary element that operates on rows, columns and data of a
1133 // given field
1134 CHKERR m_field.add_finite_element("SPRING", MF_ZERO);
1138 if (m_field.check_field(mesh_nodals_positions)) {
1140 mesh_nodals_positions);
1141 }
1142 // Add entities to that element, here we add all triangles with SPRING_BC
1143 // from cubit
1145 if (bit->getName().compare(0, 9, "SPRING_BC") == 0) {
1146 CHKERR m_field.add_ents_to_finite_element_by_dim(bit->getMeshset(), 2,
1147 "SPRING");
1148 }
1149 }
1150
1152}
1153
1155 MoFEM::Interface &m_field, const std::string field_name,
1156 const std::string mesh_nodals_positions, Range &spring_triangles) {
1158
1159 // Define boundary element that operates on rows, columns and data of a
1160 // given field
1161 CHKERR m_field.add_finite_element("SPRING_ALE", MF_ZERO);
1165 CHKERR m_field.modify_finite_element_add_field_row("SPRING_ALE",
1166 mesh_nodals_positions);
1167 CHKERR m_field.modify_finite_element_add_field_col("SPRING_ALE",
1168 mesh_nodals_positions);
1169 CHKERR m_field.modify_finite_element_add_field_data("SPRING_ALE",
1170 mesh_nodals_positions);
1171
1172 CHKERR m_field.add_ents_to_finite_element_by_dim(spring_triangles, 2,
1173 "SPRING_ALE");
1174
1176}
1177
1179 MoFEM::Interface &m_field,
1180 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr,
1181 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr,
1182 const std::string field_name, const std::string mesh_nodals_positions, double stiffness_scale) {
1184
1185 // Push operators to instances for springs
1186 // loop over blocks
1187 boost::shared_ptr<MetaSpringBC::DataAtIntegrationPtsSprings> commonDataPtr =
1188 boost::make_shared<MetaSpringBC::DataAtIntegrationPtsSprings>(m_field, stiffness_scale);
1189 CHKERR commonDataPtr->getParameters();
1190
1191 if (m_field.check_field(mesh_nodals_positions)) {
1192 CHKERR addHOOpsFace3D(mesh_nodals_positions, *fe_spring_lhs_ptr, false,
1193 false);
1194 CHKERR addHOOpsFace3D(mesh_nodals_positions, *fe_spring_rhs_ptr, false,
1195 false);
1196 }
1197
1198 for (auto &sitSpring : commonDataPtr->mapSpring) {
1199
1200 fe_spring_lhs_ptr->getOpPtrVector().push_back(
1201 new OpGetTangentSpEle(mesh_nodals_positions, commonDataPtr));
1202 fe_spring_lhs_ptr->getOpPtrVector().push_back(
1203 new OpGetNormalSpEle(mesh_nodals_positions, commonDataPtr));
1204 fe_spring_lhs_ptr->getOpPtrVector().push_back(
1205 new OpSpringKs(commonDataPtr, sitSpring.second, field_name));
1206
1207 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1208 new OpCalculateVectorFieldValues<3>(field_name, commonDataPtr->xAtPts));
1209 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1210 new OpCalculateVectorFieldValues<3>(mesh_nodals_positions,
1211 commonDataPtr->xInitAtPts));
1212 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1213 new OpGetTangentSpEle(mesh_nodals_positions, commonDataPtr));
1214 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1215 new OpGetNormalSpEle(mesh_nodals_positions, commonDataPtr));
1216 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1217 new OpSpringFs(commonDataPtr, sitSpring.second, field_name));
1218 }
1219
1221}
1222
1224 MoFEM::Interface &m_field,
1225 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr_dx,
1226 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr_dX,
1227 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr,
1228 boost::shared_ptr<MetaSpringBC::DataAtIntegrationPtsSprings>
1229 data_at_integration_pts,
1230 const std::string field_name, const std::string mesh_nodals_positions,
1231 std::string side_fe_name) {
1233
1234 // Push operators to instances for springs
1235 // loop over blocks
1236 CHKERR data_at_integration_pts->getParameters();
1237
1238 if (m_field.check_field(mesh_nodals_positions)) {
1239 CHKERR addHOOpsFace3D(mesh_nodals_positions, *fe_spring_lhs_ptr_dx, false,
1240 false);
1241 CHKERR addHOOpsFace3D(mesh_nodals_positions, *fe_spring_lhs_ptr_dX, false,
1242 false);
1243 CHKERR addHOOpsFace3D(mesh_nodals_positions, *fe_spring_rhs_ptr, false,
1244 false);
1245 }
1246
1247 for (auto &sitSpring : data_at_integration_pts->mapSpring) {
1248
1249 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> feMatSideRhs =
1250 boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(m_field);
1251
1252 feMatSideRhs->getOpPtrVector().push_back(
1254 mesh_nodals_positions, data_at_integration_pts->HMat));
1255 feMatSideRhs->getOpPtrVector().push_back(
1257 field_name, data_at_integration_pts->hMat));
1258
1259 feMatSideRhs->getOpPtrVector().push_back(new OpCalculateDeformation(
1260 mesh_nodals_positions, data_at_integration_pts));
1261
1262 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1263 new OpGetTangentSpEle(mesh_nodals_positions, data_at_integration_pts));
1264 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1265 new OpGetNormalSpEle(mesh_nodals_positions, data_at_integration_pts));
1266
1267 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1269 data_at_integration_pts->xAtPts));
1270 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1272 mesh_nodals_positions, data_at_integration_pts->xInitAtPts));
1273
1274 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1275 new OpSpringFsMaterial(mesh_nodals_positions, data_at_integration_pts,
1276 feMatSideRhs, side_fe_name, sitSpring.second));
1277
1278 fe_spring_lhs_ptr_dx->getOpPtrVector().push_back(
1280 data_at_integration_pts->xAtPts));
1281 fe_spring_lhs_ptr_dx->getOpPtrVector().push_back(
1283 mesh_nodals_positions, data_at_integration_pts->xInitAtPts));
1284 fe_spring_lhs_ptr_dx->getOpPtrVector().push_back(
1285 new OpGetTangentSpEle(mesh_nodals_positions, data_at_integration_pts));
1286 fe_spring_lhs_ptr_dx->getOpPtrVector().push_back(
1287 new OpGetNormalSpEle(mesh_nodals_positions, data_at_integration_pts));
1288
1289 fe_spring_lhs_ptr_dx->getOpPtrVector().push_back(
1290 new OpSpringKs_dX(data_at_integration_pts, sitSpring.second, field_name,
1291 mesh_nodals_positions));
1292
1293 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> feMatSideLhs_dx =
1294 boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(m_field);
1295
1296 feMatSideLhs_dx->getOpPtrVector().push_back(
1298 mesh_nodals_positions, data_at_integration_pts->HMat));
1299 feMatSideLhs_dx->getOpPtrVector().push_back(
1301 field_name, data_at_integration_pts->hMat));
1302
1303 feMatSideLhs_dx->getOpPtrVector().push_back(new OpCalculateDeformation(
1304 mesh_nodals_positions, data_at_integration_pts));
1305
1306 feMatSideLhs_dx->getOpPtrVector().push_back(
1308 mesh_nodals_positions, field_name, data_at_integration_pts));
1309
1310 fe_spring_lhs_ptr_dX->getOpPtrVector().push_back(
1312 data_at_integration_pts->xAtPts));
1313 fe_spring_lhs_ptr_dX->getOpPtrVector().push_back(
1315 mesh_nodals_positions, data_at_integration_pts->xInitAtPts));
1316 fe_spring_lhs_ptr_dX->getOpPtrVector().push_back(
1317 new OpGetTangentSpEle(mesh_nodals_positions, data_at_integration_pts));
1318 fe_spring_lhs_ptr_dX->getOpPtrVector().push_back(
1319 new OpGetNormalSpEle(mesh_nodals_positions, data_at_integration_pts));
1320
1321 fe_spring_lhs_ptr_dX->getOpPtrVector().push_back(
1322 new OpSpringALEMaterialLhs_dX_dx(mesh_nodals_positions, field_name,
1323 data_at_integration_pts,
1324 feMatSideLhs_dx, side_fe_name));
1325
1326 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> feMatSideLhs_dX =
1327 boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(m_field);
1328
1329 feMatSideLhs_dX->getOpPtrVector().push_back(
1331 mesh_nodals_positions, data_at_integration_pts->HMat));
1332 feMatSideLhs_dX->getOpPtrVector().push_back(
1334 field_name, data_at_integration_pts->hMat));
1335
1336 feMatSideLhs_dX->getOpPtrVector().push_back(new OpCalculateDeformation(
1337 mesh_nodals_positions, data_at_integration_pts));
1338
1339 feMatSideLhs_dX->getOpPtrVector().push_back(
1340 new SpringALEMaterialVolOnSideLhs_dX_dX(mesh_nodals_positions,
1341 mesh_nodals_positions,
1342 data_at_integration_pts));
1343
1344 fe_spring_lhs_ptr_dX->getOpPtrVector().push_back(
1346 mesh_nodals_positions, mesh_nodals_positions,
1347 data_at_integration_pts, feMatSideLhs_dX, side_fe_name));
1348 }
1349
1351}
Kronecker Delta class symmetric.
Kronecker Delta class.
@ 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
@ BLOCKSET
Definition: definitions.h:148
#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
FTensor::Index< 'm', SPACE_DIM > m
constexpr auto t_kd
virtual MoFEMErrorCode add_ents_to_finite_element_by_dim(const EntityHandle entities, const int dim, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual bool check_field(const std::string &name) const =0
check if field is in database
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
auto bit
set bit
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
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:1486
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1475
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
constexpr auto field_name
Operator for computing deformation gradients in side volumes.
MoFEMErrorCode doWork(int side, EntityType type, EntData &row_data)
Computes, for material configuration, normal to face that lies on a surface with springs.
MoFEMErrorCode doWork(int side, EntityType type, EntData &row_data)
Evaluates normal vector to the triangle on surface with springs based on material base coordinates.
Computes, for material configuration, tangent vectors to face that lie on a surface with springs.
MoFEMErrorCode doWork(int side, EntityType type, EntData &row_data)
Evaluates the two tangent vector to the triangle on surface with springs based on material base coord...
LHS-operator for the pressure element (material configuration)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Compute part of the left-hand side.
LHS-operator for the pressure element (material configuration)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data)
RHS-operator for the spring boundary condition element.
MetaSpringBC::BlockOptionDataSprings & dAta
boost::shared_ptr< MetaSpringBC::DataAtIntegrationPtsSprings > commonDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Integrates and assembles to global RHS vector virtual work of springs.
RHS-operator for the spring boundary condition element for ALE formulation.
MoFEMErrorCode aSsemble(EntData &row_data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode iNtegrate(EntData &row_data)
Integrates and assembles to global RHS vector virtual work of springs on material positions for confi...
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Compute part of the left-hand side.
LHS-operator for the springs element.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Integrates and assembles to global RHS vector virtual work of springs.
LHS-operator (material configuration) on the side volume.
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrates over a face contribution from a side volume.
LHS-operator (material configuration) on the side volume for spring element.
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrates over a face contribution from a side volume.
MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
static MoFEMErrorCode setSpringOperators(MoFEM::Interface &m_field, boost::shared_ptr< FaceElementForcesAndSourcesCore > fe_spring_lhs_ptr, boost::shared_ptr< FaceElementForcesAndSourcesCore > fe_spring_rhs_ptr, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS", double stiffness_scale=1.)
Implementation of spring element. Set operators to calculate LHS and RHS.
static MoFEMErrorCode addSpringElementsALE(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions, Range &spring_triangles)
Declare spring element in ALE.
static MoFEMErrorCode addSpringElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Declare spring element.
static MoFEMErrorCode setSpringOperatorsMaterial(MoFEM::Interface &m_field, boost::shared_ptr< FaceElementForcesAndSourcesCore > fe_spring_lhs_ptr_dx, boost::shared_ptr< FaceElementForcesAndSourcesCore > fe_spring_lhs_ptr_dX, boost::shared_ptr< FaceElementForcesAndSourcesCore > fe_spring_rhs_ptr, boost::shared_ptr< DataAtIntegrationPtsSprings > data_at_integration_pts, const std::string field_name, const std::string mesh_nodals_positions, std::string side_fe_name)
Implementation of spring element. Set operators to calculate LHS and RHS.
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 & 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.
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
auto getFTensor0IntegrationWeight()
Get integration weights.
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Vec & snes_f
residual