v0.13.1
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/* This file is part of MoFEM.
6 * MoFEM is free software: you can redistribute it and/or modify it under
7 * the terms of the GNU Lesser General Public License as published by the
8 * Free Software Foundation, either version 3 of the License, or (at your
9 * option) any later version.
10 *
11 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
12 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
14 * License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public
17 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
18
19#include <MoFEM.hpp>
20using namespace MoFEM;
21
22#include <SpringElement.hpp>
23using namespace boost::numeric;
24
26 EntData &data) {
27
29 // check that the faces have associated degrees of freedom
30 const int nb_dofs = data.getIndices().size();
31 if (nb_dofs == 0)
33 if (dAta.tRis.find(getFEEntityHandle()) == dAta.tRis.end()) {
35 }
36 CHKERR commonDataPtr->getBlockData(dAta);
37
38 // size of force vector associated to the entity
39 // set equal to the number of degrees of freedom of associated with the
40 // entity
41 nF.resize(nb_dofs, false);
42 nF.clear();
43
44 // get number of Gauss points
45 const int nb_gauss_pts = data.getN().size1();
46
47 // get integration weights
49
50 // FTensor indices
54 const double normal_stiffness = commonDataPtr->springStiffnessNormal;
55 const double tangent_stiffness = commonDataPtr->springStiffnessTangent;
56
57 // create a 3d vector to be used as the normal to the face with length equal
58 // to the face area
60 FTensor::Tensor2<double, 3, 3> t_normal_projection;
61 FTensor::Tensor2<double, 3, 3> t_tangent_projection;
62
63 // Extract solution at Gauss points
64 auto t_solution_at_gauss_point =
65 getFTensor1FromMat<3>(*commonDataPtr->xAtPts);
66 auto t_init_solution_at_gauss_point =
67 getFTensor1FromMat<3>(*commonDataPtr->xInitAtPts);
68 FTensor::Tensor1<double, 3> t_displacement_at_gauss_point;
69
70 auto t_normal_1 = getFTensor1FromMat<3>(commonDataPtr->normalVector);
71 // loop over all Gauss points of the face
72 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
73 t_normal(i) = t_normal_1(i);
74
75 const double normal_length = std::sqrt(t_normal(k) * t_normal(k));
76 t_normal(i) = t_normal(i) / normal_length;
77 t_normal_projection(i, j) = t_normal(i) * t_normal(j);
78 t_tangent_projection(i, j) = t_normal_projection(i, j);
79 t_tangent_projection(0, 0) -= 1;
80 t_tangent_projection(1, 1) -= 1;
81 t_tangent_projection(2, 2) -= 1;
82 // Calculate the displacement at the Gauss point
83 if (is_spatial_position) { // "SPATIAL_POSITION"
84 t_displacement_at_gauss_point(i) =
85 t_solution_at_gauss_point(i) - t_init_solution_at_gauss_point(i);
86 } else { // e.g. "DISPLACEMENT" or "U"
87 t_displacement_at_gauss_point(i) = t_solution_at_gauss_point(i);
88 }
89
90 double w = t_w * normal_length; // area was constant
91 if (getNumeredEntFiniteElementPtr()->getEntType() == MBTRI) {
92 w *= 0.5;
93 }
94
95 auto t_base_func = data.getFTensor0N(gg, 0);
96
97 // create a vector t_nf whose pointer points an array of 3 pointers
98 // pointing to nF memory location of components
100 &nF[2]);
101
102 for (int rr = 0; rr != nb_dofs / 3; ++rr) { // loop over the nodes
103
104 t_nf(i) += w * t_base_func * normal_stiffness *
105 t_normal_projection(i, j) * t_displacement_at_gauss_point(j);
106 t_nf(i) -= w * t_base_func * tangent_stiffness *
107 t_tangent_projection(i, j) * t_displacement_at_gauss_point(j);
108
109 // move to next base function
110 ++t_base_func;
111 // move the pointer to next element of t_nf
112 ++t_nf;
113 }
114 // move to next integration weight
115 ++t_w;
116 // move to the solutions at the next Gauss point
117 ++t_solution_at_gauss_point;
118 ++t_init_solution_at_gauss_point;
119 ++t_normal_1;
120 }
121 // add computed values of spring in the global right hand side vector
122 Vec f = getFEMethod()->ksp_f != PETSC_NULL ? getFEMethod()->ksp_f
123 : getFEMethod()->snes_f;
124 CHKERR VecSetValues(f, nb_dofs, &data.getIndices()[0], &nF[0], ADD_VALUES);
126}
127
129 EntityType row_type,
130 EntityType col_type,
131 EntData &row_data,
132 EntData &col_data) {
134 // check if the volumes have associated degrees of freedom
135 const int row_nb_dofs = row_data.getIndices().size();
136 if (!row_nb_dofs)
138
139 const int col_nb_dofs = col_data.getIndices().size();
140 if (!col_nb_dofs)
142
143 if (dAta.tRis.find(getFEEntityHandle()) == dAta.tRis.end()) {
145 }
146
147 CHKERR commonDataPtr->getBlockData(dAta);
148 // size associated to the entity
149 locKs.resize(row_nb_dofs, col_nb_dofs, false);
150 locKs.clear();
151
152 // get number of Gauss points
153 const int row_nb_gauss_pts = row_data.getN().size1();
154 if (!row_nb_gauss_pts) // check if number of Gauss point <> 0
156
157 // get intergration weights
158 auto t_w = getFTensor0IntegrationWeight();
159
163
164 // create a 3d vector to be used as the normal to the face with length equal
165 // to the face area
167
168 // First tangent vector
169 auto t_tangent1_ptr = getFTensor1Tangent1AtGaussPts();
171 t_tangent1(i) = t_tangent1_ptr(i);
172
173 // Second tangent vector, such that t_n = t_t1 x t_t2 | t_t2 = t_n x t_t1
175 t_tangent2(i) = FTensor::levi_civita(i, j, k) * t_normal(j) * t_tangent1(k);
176
177 auto t_tangent2_ptr = getFTensor1Tangent2AtGaussPts();
178 t_normal(i) =
179 FTensor::levi_civita(i, j, k) * t_tangent1(j) * t_tangent2_ptr(k);
180
181 auto normal_original_slave =
182 getFTensor1FromMat<3>(commonDataPtr->normalVector);
183
184 FTensor::Tensor2<double, 3, 3> t_normal_projection;
185 FTensor::Tensor2<double, 3, 3> t_tangent_projection;
186 // loop over the Gauss points
187 for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
188 t_normal(i) = normal_original_slave(i);
189 const double normal_length = std::sqrt(t_normal(k) * t_normal(k));
190 t_normal(i) = t_normal(i) / normal_length;
191 t_normal_projection(i, j) = t_normal(i) * t_normal(j);
192 t_tangent_projection(i, j) = t_normal_projection(i, j);
193 t_tangent_projection(0, 0) -= 1;
194 t_tangent_projection(1, 1) -= 1;
195 t_tangent_projection(2, 2) -= 1;
196
197 // get area and integration weight
198 double w = t_w * normal_length;
199 if (getNumeredEntFiniteElementPtr()->getEntType() == MBTRI) {
200 w *= 0.5;
201 }
202
203 auto t_row_base_func = row_data.getFTensor0N(gg, 0);
204
205 for (int rr = 0; rr != row_nb_dofs / 3; rr++) {
206 auto assemble_m = getFTensor2FromArray<3, 3, 3>(locKs, 3 * rr);
207 auto t_col_base_func = col_data.getFTensor0N(gg, 0);
208 for (int cc = 0; cc != col_nb_dofs / 3; cc++) {
209 assemble_m(i, j) += w * t_row_base_func * t_col_base_func *
210 commonDataPtr->springStiffnessNormal *
211 t_normal_projection(i, j);
212 assemble_m(i, j) -= w * t_row_base_func * t_col_base_func *
213 commonDataPtr->springStiffnessTangent *
214 t_tangent_projection(i, j);
215 ++t_col_base_func;
216 ++assemble_m;
217 }
218 ++t_row_base_func;
219 }
220 // move to next integration weight
221 ++t_w;
222 ++normal_original_slave;
223 }
224
225 // Add computed values of spring stiffness to the global LHS matrix
226 CHKERR MatSetValues(getKSPB(), row_data, col_data, &locKs(0, 0), ADD_VALUES);
227
228 // is symmetric
229 if (row_side != col_side || row_type != col_type) {
230 transLocKs.resize(col_nb_dofs, row_nb_dofs, false);
231 noalias(transLocKs) = trans(locKs);
232
233 CHKERR MatSetValues(getKSPB(), col_data, row_data, &transLocKs(0, 0),
234 ADD_VALUES);
235 }
236
238}
239
241 EntityType row_type,
242 EntityType col_type,
243 EntData &row_data,
244 EntData &col_data) {
246
247 // check if the volumes have associated degrees of freedom
248 const int row_nb_dofs = row_data.getIndices().size();
249 if (!row_nb_dofs)
251
252 const int col_nb_dofs = col_data.getIndices().size();
253 if (!col_nb_dofs)
255
256 if (dAta.tRis.find(getFEEntityHandle()) == dAta.tRis.end()) {
258 }
259
260 CHKERR commonDataPtr->getBlockData(dAta);
261 // size associated to the entity
262 locKs.resize(row_nb_dofs, col_nb_dofs, false);
263 locKs.clear();
264
265 // get number of Gauss points
266 const int row_nb_gauss_pts = row_data.getN().size1();
267 if (!row_nb_gauss_pts) // check if number of Gauss point <> 0
269
270 // get intergration weights
271 auto t_w = getFTensor0IntegrationWeight();
272
277
278 auto make_vec_der = [&](auto t_N, auto t_1, auto t_2) {
280 t_n(i, j) = 0;
281 t_n(i, j) += FTensor::levi_civita(i, j, k) * t_2(k) * t_N(0);
282 t_n(i, j) -= FTensor::levi_civita(i, j, k) * t_1(k) * t_N(1);
283 return t_n;
284 };
285
286 auto make_vec_der_2 = [&](auto t_N, auto t_1, auto t_2) {
288 t_normal(i) = FTensor::levi_civita(i, j, k) * t_1(j) * t_2(k);
289 const double normal_norm = std::sqrt(t_normal(i) * t_normal(i));
291 t_n(i, j) = 0;
292 t_n(i, j) += FTensor::levi_civita(i, j, k) * t_2(k) * t_N(0);
293 t_n(i, j) -= FTensor::levi_civita(i, j, k) * t_1(k) * t_N(1);
295 t_result(i, j) = 0;
296 t_result(i, j) =
297 t_n(i, j) / normal_norm - t_normal(i) * t_n(k, j) * t_normal(k) /
298 (normal_norm * normal_norm * normal_norm);
299 return t_result;
300 };
301
302 const double normal_stiffness = commonDataPtr->springStiffnessNormal;
303 const double tangent_stiffness = commonDataPtr->springStiffnessTangent;
304
305 // create a 3d vector to be used as the normal to the face with length equal
306 // to the face area
308 FTensor::Tensor2<double, 3, 3> t_normal_projection;
309 FTensor::Tensor2<double, 3, 3> t_tangent_projection;
310 FTensor::Tensor1<double, 3> t_displacement_at_gauss_point;
311
312 // Extract solution at Gauss points
313 auto t_solution_at_gauss_point =
314 getFTensor1FromMat<3>(*commonDataPtr->xAtPts);
315 auto t_init_solution_at_gauss_point =
316 getFTensor1FromMat<3>(*commonDataPtr->xInitAtPts);
317
318 auto t_1 = getFTensor1FromMat<3>(commonDataPtr->tangent1);
319 auto t_2 = getFTensor1FromMat<3>(commonDataPtr->tangent2);
320
321 constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
322
323 auto normal_at_gp = getFTensor1FromMat<3>(commonDataPtr->normalVector);
324
325 // loop over the Gauss points
326 for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
327 // get area and integration weight
328 t_normal(i) = normal_at_gp(i);
329 const double normal_length = std::sqrt(t_normal(k) * t_normal(k));
330 t_normal(i) = t_normal(i) / normal_length;
331
332 t_normal_projection(i, j) = t_normal(i) * t_normal(j);
333 t_tangent_projection(i, j) = t_normal_projection(i, j);
334 t_tangent_projection(0, 0) -= 1;
335 t_tangent_projection(1, 1) -= 1;
336 t_tangent_projection(2, 2) -= 1;
337
338 double w = 0.5 * t_w;
339
340 t_displacement_at_gauss_point(i) =
341 t_solution_at_gauss_point(i) - t_init_solution_at_gauss_point(i);
342
343 auto t_row_base_func = row_data.getFTensor0N(gg, 0);
344
345 for (int rr = 0; rr != row_nb_dofs / 3; rr++) {
346 auto t_col_base_func = col_data.getFTensor0N(gg, 0);
347 auto t_N = col_data.getFTensor1DiffN<2>(gg, 0);
348 auto assemble_m = getFTensor2FromArray<3, 3, 3>(locKs, 3 * rr);
349 for (int cc = 0; cc != col_nb_dofs / 3; cc++) {
350 auto t_d_n = make_vec_der(t_N, t_1, t_2);
351 auto t_d_n_2 = make_vec_der_2(t_N, t_1, t_2);
352
353 assemble_m(i, j) -= w * normal_length * t_col_base_func *
354 normal_stiffness * t_row_base_func *
355 t_normal_projection(i, j);
356
357 assemble_m(i, j) += w * normal_length * t_col_base_func *
358 tangent_stiffness * t_row_base_func *
359 t_tangent_projection(i, j);
360
361 assemble_m(i, j) +=
362 w * (normal_stiffness - tangent_stiffness) * t_row_base_func *
363 (t_d_n(i, j) * (t_normal(k) * t_displacement_at_gauss_point(k)) +
364 normal_length * t_normal(i) *
365 (t_d_n_2(k, j) * t_displacement_at_gauss_point(k)));
366 // constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
367 assemble_m(i, j) += w * (t_d_n(l, j) * t_normal(l)) *
368 tangent_stiffness * t_row_base_func * t_kd(i, k) *
369 t_displacement_at_gauss_point(k);
370
371 ++t_col_base_func;
372 ++t_N;
373 ++assemble_m;
374 }
375 ++t_row_base_func;
376 }
377 // move to next integration weight
378 ++t_w;
379 ++t_solution_at_gauss_point;
380 ++t_init_solution_at_gauss_point;
381 ++t_1;
382 ++t_2;
383 ++normal_at_gp;
384 }
385
386 // Add computed values of spring stiffness to the global LHS matrix
387 CHKERR MatSetValues(getKSPB(), row_data, col_data, &locKs(0, 0), ADD_VALUES);
388
390}
391
393 int row_side, int col_side, EntityType row_type, EntityType col_type,
394 EntData &row_data, EntData &col_data) {
395
397
398 if (dataAtSpringPts->faceRowData == nullptr)
400
401 if (row_type != MBVERTEX)
403
404 row_nb_dofs = dataAtSpringPts->faceRowData->getIndices().size();
405 if (!row_nb_dofs)
407 col_nb_dofs = col_data.getIndices().size();
408 if (!col_nb_dofs)
410
411 nb_gauss_pts = dataAtSpringPts->faceRowData->getN().size1();
412
413 nb_base_fun_row = dataAtSpringPts->faceRowData->getFieldData().size() / 3;
414 nb_base_fun_col = col_data.getFieldData().size() / 3;
415
416 NN.resize(3 * nb_base_fun_row, 3 * nb_base_fun_col, false);
417 NN.clear();
418
419 // integrate local matrix for entity block
420 CHKERR iNtegrate(*(dataAtSpringPts->faceRowData), col_data);
421
422 // assemble local matrix
423 CHKERR aSsemble(*(dataAtSpringPts->faceRowData), col_data);
424
426}
427
430 EntData &col_data) {
431
433
434 // get pointer to first global index on row
435 const int *row_indices = &*row_data.getIndices().data().begin();
436 // get pointer to first global index on column
437 const int *col_indices = &*col_data.getIndices().data().begin();
438
439 auto &data = *dataAtSpringPts;
440 if (!data.forcesOnlyOnEntitiesRow.empty()) {
441 rowIndices.resize(row_nb_dofs, false);
442 noalias(rowIndices) = row_data.getIndices();
443 row_indices = &rowIndices[0];
444 VectorDofs &dofs = row_data.getFieldDofs();
445 VectorDofs::iterator dit = dofs.begin();
446 for (int ii = 0; dit != dofs.end(); ++dit, ++ii) {
447 if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
448 data.forcesOnlyOnEntitiesRow.end()) {
449 rowIndices[ii] = -1;
450 }
451 }
452 }
453
454 if (!data.forcesOnlyOnEntitiesCol.empty()) {
455 colIndices.resize(col_nb_dofs, false);
456 noalias(colIndices) = col_data.getIndices();
457 col_indices = &colIndices[0];
458 VectorDofs &dofs = col_data.getFieldDofs();
459 VectorDofs::iterator dit = dofs.begin();
460 for (int ii = 0; dit != dofs.end(); ++dit, ++ii) {
461 if (data.forcesOnlyOnEntitiesCol.find((*dit)->getEnt()) ==
462 data.forcesOnlyOnEntitiesCol.end()) {
463 colIndices[ii] = -1;
464 }
465 }
466 }
467
468 Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
469 : getFEMethod()->snes_B;
470 // assemble local matrix
471 CHKERR MatSetValues(B, row_nb_dofs, row_indices, col_nb_dofs, col_indices,
472 &*NN.data().begin(), ADD_VALUES);
473
475}
476
478 EntData &row_data, EntData &col_data) {
479
481
486
487 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
489 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
490 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
491 &m(r + 2, c + 2));
492 };
494 FTensor::Tensor2<double, 3, 3> t_normal_projection;
495 FTensor::Tensor2<double, 3, 3> t_tangent_projection;
496
497 const double normal_stiffness = dataAtSpringPts->springStiffnessNormal;
498 const double tangent_stiffness = dataAtSpringPts->springStiffnessTangent;
499
500 // Extract solution at Gauss points
501 auto t_solution_at_gauss_point =
502 getFTensor1FromMat<3>(*dataAtSpringPts->xAtPts);
503 auto t_init_solution_at_gauss_point =
504 getFTensor1FromMat<3>(*dataAtSpringPts->xInitAtPts);
505 FTensor::Tensor1<double, 3> t_displacement_at_gauss_point;
506
507 auto t_normal_1 = getFTensor1FromMat<3>(dataAtSpringPts->normalVector);
508 auto t_w = getFTensor0IntegrationWeight();
509
510 auto t_inv_H = getFTensor2FromMat<3, 3>(*dataAtSpringPts->invHMat);
511 auto t_F = getFTensor2FromMat<3, 3>(*dataAtSpringPts->FMat);
512 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
513
514 t_normal(i) = t_normal_1(i);
515
516 const double normal_length = std::sqrt(t_normal(k) * t_normal(k));
517 t_normal(i) = t_normal(i) / normal_length;
518 t_normal_projection(i, j) = t_normal(i) * t_normal(j);
519 t_tangent_projection(i, j) = t_normal_projection(i, j);
520 t_tangent_projection(0, 0) -= 1;
521 t_tangent_projection(1, 1) -= 1;
522 t_tangent_projection(2, 2) -= 1;
523
524 t_displacement_at_gauss_point(i) =
525 t_solution_at_gauss_point(i) - t_init_solution_at_gauss_point(i);
526
527 double w = 0.5 * t_w * normal_length;
528
529 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
530 FTensor::Tensor0<double *> t_col_base(&col_data.getN()(gg, 0));
531 int bbc = 0;
532 for (; bbc != nb_base_fun_col; ++bbc) {
533
534 FTensor::Tensor0<double *> t_row_base(&row_data.getN()(gg, 0));
535
536 int bbr = 0;
537 for (; bbr != nb_base_fun_row; ++bbr) {
538
539 auto t_assemble = get_tensor2(NN, 3 * bbr, 3 * bbc);
540
541 t_assemble(i, j) -= w * t_row_base * t_col_base * normal_stiffness *
542 t_F(k, i) * t_normal_projection(k, j);
543 t_assemble(i, j) += w * t_row_base * t_col_base * tangent_stiffness *
544 t_F(k, i) * t_tangent_projection(k, j);
545
546 t_assemble(i, j) -= w * normal_stiffness * t_row_base * t_inv_H(k, i) *
547 t_col_diff_base(k) * t_normal_projection(j, l) *
548 t_displacement_at_gauss_point(l);
549 t_assemble(i, j) += w * tangent_stiffness * t_row_base * t_inv_H(k, i) *
550 t_col_diff_base(k) * t_tangent_projection(j, l) *
551 t_displacement_at_gauss_point(l);
552
553 ++t_row_base;
554 }
555 ++t_col_diff_base;
556 ++t_col_base;
557 }
558 ++t_w;
559 ++t_solution_at_gauss_point;
560 ++t_init_solution_at_gauss_point;
561 ++t_normal_1;
562 ++t_inv_H;
563 ++t_F;
564 }
565
567}
568
571 EntData &col_data) {
572
574
575 // get pointer to first global index on row
576 const int *row_indices = &*row_data.getIndices().data().begin();
577 // get pointer to first global index on column
578 const int *col_indices = &*col_data.getIndices().data().begin();
579
580 auto &data = *dataAtSpringPts;
581 if (!data.forcesOnlyOnEntitiesRow.empty()) {
582 rowIndices.resize(row_nb_dofs, false);
583 noalias(rowIndices) = row_data.getIndices();
584 row_indices = &rowIndices[0];
585 VectorDofs &dofs = row_data.getFieldDofs();
586 VectorDofs::iterator dit = dofs.begin();
587 for (int ii = 0; dit != dofs.end(); ++dit, ++ii) {
588 if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
589 data.forcesOnlyOnEntitiesRow.end()) {
590 rowIndices[ii] = -1;
591 }
592 }
593 }
594
595 if (!data.forcesOnlyOnEntitiesCol.empty()) {
596 colIndices.resize(col_nb_dofs, false);
597 noalias(colIndices) = col_data.getIndices();
598 col_indices = &colIndices[0];
599 VectorDofs &dofs = col_data.getFieldDofs();
600 VectorDofs::iterator dit = dofs.begin();
601 for (int ii = 0; dit != dofs.end(); ++dit, ++ii) {
602 if (data.forcesOnlyOnEntitiesCol.find((*dit)->getEnt()) ==
603 data.forcesOnlyOnEntitiesCol.end()) {
604 colIndices[ii] = -1;
605 }
606 }
607 }
608
609 Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
610 : getFEMethod()->snes_B;
611 // assemble local matrix
612 CHKERR MatSetValues(B, row_nb_dofs, row_indices, col_nb_dofs, col_indices,
613 &*NN.data().begin(), ADD_VALUES);
614
616}
617
619 int row_side, int col_side, EntityType row_type, EntityType col_type,
620 EntData &row_data, EntData &col_data) {
621
623
624 row_nb_dofs = row_data.getIndices().size();
625 if (!row_nb_dofs)
627 col_nb_dofs = col_data.getIndices().size();
628 if (!col_nb_dofs)
630 nb_gauss_pts = row_data.getN().size1();
631
632 nb_base_fun_row = row_data.getFieldData().size() / 3;
633 nb_base_fun_col = col_data.getFieldData().size() / 3;
634
635 NN.resize(3 * nb_base_fun_row, 3 * nb_base_fun_col, false);
636 NN.clear();
637
638 diagonal_block = (row_type == col_type) && (row_side == col_side);
639
640 if (col_type == MBVERTEX) {
641 dataAtSpringPts->faceRowData = &row_data;
642 CHKERR loopSideVolumes(sideFeName, *sideFe);
643 }
644
645 // integrate local matrix for entity block
646 CHKERR iNtegrate(row_data, col_data);
647
648 // assemble local matrix
649 CHKERR aSsemble(row_data, col_data);
650
652}
653
656 EntData &col_data) {
657
659
664
665 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
667 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
668 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
669 &m(r + 2, c + 2));
670 };
671
672 auto make_vec_der = [&](auto t_N, auto t_1, auto t_2) {
674 t_n(i, j) = 0;
675 t_n(i, j) += FTensor::levi_civita(i, j, k) * t_2(k) * t_N(0);
676 t_n(i, j) -= FTensor::levi_civita(i, j, k) * t_1(k) * t_N(1);
677 return t_n;
678 };
679
680 auto make_vec_der_2 = [&](auto t_N, auto t_1, auto t_2) {
682 t_normal(i) = FTensor::levi_civita(i, j, k) * t_1(j) * t_2(k);
683 const double normal_norm = std::sqrt(t_normal(i) * t_normal(i));
685 t_n(i, j) = 0;
686 t_n(i, j) += FTensor::levi_civita(i, j, k) * t_2(k) * t_N(0);
687 t_n(i, j) -= FTensor::levi_civita(i, j, k) * t_1(k) * t_N(1);
689 t_result(i, j) = 0;
690 t_result(i, j) =
691 t_n(i, j) / normal_norm - t_normal(i) * t_n(k, j) * t_normal(k) /
692 (normal_norm * normal_norm * normal_norm);
693 return t_result;
694 };
695
696 dataAtSpringPts->faceRowData = nullptr;
697 CHKERR loopSideVolumes(sideFeName, *sideFe);
698
699 auto t_F = getFTensor2FromMat<3, 3>(*dataAtSpringPts->FMat);
700
701 auto t_w = getFTensor0IntegrationWeight();
702 auto t_1 = getFTensor1FromMat<3>(dataAtSpringPts->tangent1);
703 auto t_2 = getFTensor1FromMat<3>(dataAtSpringPts->tangent2);
704
705 const double normal_stiffness = dataAtSpringPts->springStiffnessNormal;
706 const double tangent_stiffness = dataAtSpringPts->springStiffnessTangent;
707
708 // create a 3d vector to be used as the normal to the face with length equal
709 // to the face area
711 FTensor::Tensor2<double, 3, 3> t_normal_projection;
712 FTensor::Tensor2<double, 3, 3> t_tangent_projection;
713
714 // Extract solution at Gauss points
715 auto t_solution_at_gauss_point =
716 getFTensor1FromMat<3>(*dataAtSpringPts->xAtPts);
717 auto t_init_solution_at_gauss_point =
718 getFTensor1FromMat<3>(*dataAtSpringPts->xInitAtPts);
719 FTensor::Tensor1<double, 3> t_displacement_at_gauss_point;
720
721 auto t_normal_1 = getFTensor1FromMat<3>(dataAtSpringPts->normalVector);
723
724 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
725 t_normal(i) = t_normal_1(i);
726
727 const double normal_length = std::sqrt(t_normal(k) * t_normal(k));
728 t_normal(i) = t_normal(i) / normal_length;
729 t_normal_projection(i, j) = t_normal(i) * t_normal(j);
730 t_tangent_projection(i, j) = t_normal_projection(i, j);
731 t_tangent_projection(0, 0) -= 1;
732 t_tangent_projection(1, 1) -= 1;
733 t_tangent_projection(2, 2) -= 1;
734 // Calculate the displacement at the Gauss point
735 t_displacement_at_gauss_point(i) =
736 t_solution_at_gauss_point(i) - t_init_solution_at_gauss_point(i);
737 double val = 0.5 * t_w;
738
739 auto t_N = col_data.getFTensor1DiffN<2>(gg, 0);
740 auto t_col_base_func = col_data.getFTensor0N(gg, 0);
741 int bbc = 0;
742 for (; bbc != nb_base_fun_col; ++bbc) {
743
744 auto t_row_base_func = row_data.getFTensor0N(gg, 0);
745
746 int bbr = 0;
747 for (; bbr != nb_base_fun_row; ++bbr) {
748
749 auto t_d_n = make_vec_der(t_N, t_1, t_2);
750
751 auto t_assemble = get_tensor2(NN, 3 * bbr, 3 * bbc);
752
753 auto t_d_n_2 = make_vec_der_2(t_N, t_1, t_2);
754
755 t_assemble(i, j) += val * normal_length * t_col_base_func *
756 normal_stiffness * t_row_base_func * t_F(k, i) *
757 t_normal_projection(k, j);
758
759 t_assemble(i, j) -= val * normal_length * t_col_base_func *
760 tangent_stiffness * t_row_base_func * t_F(k, i) *
761 t_tangent_projection(k, j);
762
763 t_assemble(i, j) -=
764 val * (normal_stiffness - tangent_stiffness) * t_row_base_func *
765 t_F(l, i) *
766 (t_d_n(l, j) * (t_normal(k) * t_displacement_at_gauss_point(k)) +
767 normal_length * t_normal(l) *
768 (t_d_n_2(k, j) * t_displacement_at_gauss_point(k)));
769 // constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
770 t_assemble(i, j) -= val * (t_d_n(l, j) * t_normal(l)) *
771 tangent_stiffness * t_row_base_func * t_F(k, i) *
772 t_displacement_at_gauss_point(k);
773
774 ++t_row_base_func;
775 }
776 ++t_N;
777 ++t_col_base_func;
778 }
779 ++t_F;
780 ++t_w;
781 ++t_1;
782 ++t_2;
783 ++t_normal_1;
784 ++t_solution_at_gauss_point;
785 ++t_init_solution_at_gauss_point;
786 }
787
789}
790
792 EntData &row_data, EntData &col_data) {
793
795
802
803 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
805 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
806 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
807 &m(r + 2, c + 2));
808 };
809
810 const double normal_stiffness = dataAtSpringPts->springStiffnessNormal;
811 const double tangent_stiffness = dataAtSpringPts->springStiffnessTangent;
812
813 // create a 3d vector to be used as the normal to the face with length equal
814 // to the face area
816 FTensor::Tensor2<double, 3, 3> t_normal_projection;
817 FTensor::Tensor2<double, 3, 3> t_tangent_projection;
818
819 // Extract solution at Gauss points
820 auto t_solution_at_gauss_point =
821 getFTensor1FromMat<3>(*dataAtSpringPts->xAtPts);
822 auto t_init_solution_at_gauss_point =
823 getFTensor1FromMat<3>(*dataAtSpringPts->xInitAtPts);
824 FTensor::Tensor1<double, 3> t_displacement_at_gauss_point;
825
826 auto t_normal_1 = getFTensor1FromMat<3>(dataAtSpringPts->normalVector);
827
828 auto t_w = getFTensor0IntegrationWeight();
829
830 auto t_h = getFTensor2FromMat<3, 3>(*dataAtSpringPts->hMat);
831 auto t_inv_H = getFTensor2FromMat<3, 3>(*dataAtSpringPts->invHMat);
832
834
835 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
836 t_normal(i) = t_normal_1(i);
837
838 const double normal_length = std::sqrt(t_normal(k) * t_normal(k));
839 t_normal(i) = t_normal(i) / normal_length;
840 t_normal_projection(i, j) = t_normal(i) * t_normal(j);
841 t_tangent_projection(i, j) = t_normal_projection(i, j);
842 t_tangent_projection(0, 0) -= 1;
843 t_tangent_projection(1, 1) -= 1;
844 t_tangent_projection(2, 2) -= 1;
845 // Calculate the displacement at the Gauss point
846 t_displacement_at_gauss_point(i) =
847 t_solution_at_gauss_point(i) - t_init_solution_at_gauss_point(i);
848
849 double val = 0.5 * t_w * normal_length;
850
851 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
852
853 int bbc = 0;
854 for (; bbc != nb_base_fun_col; ++bbc) {
855
856 FTensor::Tensor0<double *> t_row_base(&row_data.getN()(gg, 0));
857
858 int bbr = 0;
859 for (; bbr != nb_base_fun_row; ++bbr) {
860
861 auto t_assemble = get_tensor2(NN, 3 * bbr, 3 * bbc);
862
863 t_assemble(i, j) +=
864 val * t_row_base * normal_stiffness * t_inv_H(l, j) *
865 t_col_diff_base(m) * t_inv_H(m, i) * t_h(k, l) *
866 (t_normal_projection(k, q) * t_displacement_at_gauss_point(q));
867
868 t_assemble(i, j) -=
869 val * t_row_base * tangent_stiffness * t_inv_H(l, j) *
870 t_col_diff_base(m) * t_inv_H(m, i) * t_h(k, l) *
871 (t_tangent_projection(k, q) * t_displacement_at_gauss_point(q));
872
873 ++t_row_base;
874 }
875 ++t_col_diff_base;
876 }
877 ++t_w;
878 ++t_h;
879 ++t_inv_H;
880 ++t_solution_at_gauss_point;
881 ++t_init_solution_at_gauss_point;
882 ++t_normal_1;
883 }
884
886}
887
890 EntData &row_data) {
891
896
897 // get number of integration points
898 const int nb_integration_pts = getGaussPts().size2();
899
900 auto t_h = getFTensor2FromMat<3, 3>(*commonDataPtr->hMat);
901 auto t_H = getFTensor2FromMat<3, 3>(*commonDataPtr->HMat);
902
903 commonDataPtr->detHVec->resize(nb_integration_pts, false);
904 commonDataPtr->invHMat->resize(9, nb_integration_pts, false);
905 commonDataPtr->FMat->resize(9, nb_integration_pts, false);
906
907 commonDataPtr->detHVec->clear();
908 commonDataPtr->invHMat->clear();
909 commonDataPtr->FMat->clear();
910
911 auto t_detH = getFTensor0FromVec(*commonDataPtr->detHVec);
912 auto t_invH = getFTensor2FromMat<3, 3>(*commonDataPtr->invHMat);
913 auto t_F = getFTensor2FromMat<3, 3>(*commonDataPtr->FMat);
914
915 for (int gg = 0; gg != nb_integration_pts; ++gg) {
916 CHKERR determinantTensor3by3(t_H, t_detH);
917 CHKERR invertTensor3by3(t_H, t_detH, t_invH);
918 t_F(i, j) = t_h(i, k) * t_invH(k, j);
919 ++t_h;
920 ++t_H;
921 ++t_detH;
922 ++t_invH;
923 ++t_F;
924 }
925
927}
928
931 EntData &row_data) {
932
934
935 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
936 dAta.tRis.end()) {
938 }
939
940 // get number of dofs on row
941 nbRows = row_data.getIndices().size();
942 // if no dofs on row, exit that work, nothing to do here
943 if (!nbRows)
945
946 nF.resize(nbRows, false);
947 nF.clear();
948
949 // get number of integration points
950 nbIntegrationPts = getGaussPts().size2();
951
952 // integrate local matrix for entity block
953 CHKERR iNtegrate(row_data);
954
955 // assemble local matrix
956 CHKERR aSsemble(row_data);
957
959}
960
963
964 CHKERR loopSideVolumes(sideFeName, *sideFe);
965
966 // check that the faces have associated degrees of freedom
967 const int nb_dofs = data.getIndices().size();
968 if (nb_dofs == 0)
970
971 CHKERR dataAtPts->getBlockData(dAta);
972
973 // get integration weights
974 auto t_w = getFTensor0IntegrationWeight();
975
976 // FTensor indices
980
981 const double normal_stiffness = dataAtPts->springStiffnessNormal;
982 const double tangent_stiffness = dataAtPts->springStiffnessTangent;
983
984 // create a 3d vector to be used as the normal to the face with length equal
985 // to the face area
987 FTensor::Tensor2<double, 3, 3> t_normal_projection;
988 FTensor::Tensor2<double, 3, 3> t_tangent_projection;
989
990 // Extract solution at Gauss points
991 auto t_solution_at_gauss_point = getFTensor1FromMat<3>(*dataAtPts->xAtPts);
992 auto t_init_solution_at_gauss_point =
993 getFTensor1FromMat<3>(*dataAtPts->xInitAtPts);
994 FTensor::Tensor1<double, 3> t_displacement_at_gauss_point;
995
996 auto t_normal_1 = getFTensor1FromMat<3>(dataAtPts->normalVector);
997
998 auto t_F = getFTensor2FromMat<3, 3>(*dataAtPts->FMat);
999
1000 // loop over all Gauss points of the face
1001 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
1002 t_normal(i) = t_normal_1(i);
1003
1004 const double normal_length = std::sqrt(t_normal(k) * t_normal(k));
1005 t_normal(i) = t_normal(i) / normal_length;
1006 t_normal_projection(i, j) = t_normal(i) * t_normal(j);
1007 t_tangent_projection(i, j) = t_normal_projection(i, j);
1008 t_tangent_projection(0, 0) -= 1;
1009 t_tangent_projection(1, 1) -= 1;
1010 t_tangent_projection(2, 2) -= 1;
1011 // Calculate the displacement at the Gauss point
1012 t_displacement_at_gauss_point(i) =
1013 t_solution_at_gauss_point(i) - t_init_solution_at_gauss_point(i);
1014
1015 double w = t_w * 0.5 * normal_length; // area was constant
1016
1017 auto t_base_func = data.getFTensor0N(gg, 0);
1018
1019 // create a vector t_nf whose pointer points an array of 3 pointers
1020 // pointing to nF memory location of components
1021 FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_nf(&nF[0], &nF[1],
1022 &nF[2]);
1023
1024 for (int rr = 0; rr != nb_dofs / 3; ++rr) { // loop over the nodes
1025
1026 t_nf(i) -= w * t_base_func * normal_stiffness * t_F(k, i) *
1027 t_normal_projection(k, j) * t_displacement_at_gauss_point(j);
1028 t_nf(i) += w * t_base_func * tangent_stiffness * t_F(k, i) *
1029 t_tangent_projection(k, j) * t_displacement_at_gauss_point(j);
1030
1031 // move to next base function
1032 ++t_base_func;
1033 // move the pointer to next element of t_nf
1034 ++t_nf;
1035 }
1036 // move to next integration weight
1037 ++t_w;
1038 // move to the solutions at the next Gauss point
1039 ++t_solution_at_gauss_point;
1040 ++t_init_solution_at_gauss_point;
1041 ++t_normal_1;
1042 ++t_F;
1043 }
1044 // add computed values of spring in the global right hand side vector
1045 // Vec f = getFEMethod()->ksp_f != PETSC_NULL ? getFEMethod()->ksp_f
1046 // : getFEMethod()->snes_f;
1047 // CHKERR VecSetValues(f, nb_dofs, &data.getIndices()[0], &nF[0], ADD_VALUES);
1049}
1050
1053
1054 // get pointer to first global index on row
1055 const int *row_indices = &*row_data.getIndices().data().begin();
1056
1057 auto &data = *dataAtPts;
1058 if (!data.forcesOnlyOnEntitiesRow.empty()) {
1059 rowIndices.resize(nbRows, false);
1060 noalias(rowIndices) = row_data.getIndices();
1061 row_indices = &rowIndices[0];
1062 VectorDofs &dofs = row_data.getFieldDofs();
1063 VectorDofs::iterator dit = dofs.begin();
1064 for (int ii = 0; dit != dofs.end(); ++dit, ++ii) {
1065 if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
1066 data.forcesOnlyOnEntitiesRow.end()) {
1067 rowIndices[ii] = -1;
1068 }
1069 }
1070 }
1071
1072 CHKERR VecSetOption(getSNESf(), VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
1073 CHKERR VecSetValues(getSNESf(), nbRows, row_indices, &*nF.data().begin(),
1074 ADD_VALUES);
1075
1077}
1078
1081 EntData &data) {
1083
1084 if (data.getFieldData().size() == 0)
1085 PetscFunctionReturn(0);
1086
1087 ngp = data.getN().size1();
1088
1089 unsigned int nb_dofs = data.getFieldData().size() / 3;
1091 if (type == MBVERTEX) {
1092 dataAtIntegrationPts->tangent1.resize(3, ngp, false);
1093 dataAtIntegrationPts->tangent1.clear();
1094
1095 dataAtIntegrationPts->tangent2.resize(3, ngp, false);
1096 dataAtIntegrationPts->tangent2.clear();
1097 }
1098
1099 auto t_1 = getFTensor1FromMat<3>(dataAtIntegrationPts->tangent1);
1100 auto t_2 = getFTensor1FromMat<3>(dataAtIntegrationPts->tangent2);
1101
1102 for (unsigned int gg = 0; gg != ngp; ++gg) {
1103 auto t_N = data.getFTensor1DiffN<2>(gg, 0);
1104 auto t_dof = data.getFTensor1FieldData<3>();
1105
1106 for (unsigned int dd = 0; dd != nb_dofs; ++dd) {
1107 t_1(i) += t_dof(i) * t_N(0);
1108 t_2(i) += t_dof(i) * t_N(1);
1109 ++t_dof;
1110 ++t_N;
1111 }
1112 ++t_1;
1113 ++t_2;
1114 }
1115
1117}
1118
1120 EntData &data) {
1122
1123 if (data.getFieldData().size() == 0)
1125
1126 ngp = data.getN().size1();
1127
1131
1132 if (type == MBVERTEX) {
1133 dataAtIntegrationPts->normalVector.resize(3, ngp, false);
1134 dataAtIntegrationPts->normalVector.clear();
1135 }
1136
1137 auto normal_original_slave =
1138 getFTensor1FromMat<3>(dataAtIntegrationPts->normalVector);
1139
1140 auto tangent_1 = getFTensor1FromMat<3>(dataAtIntegrationPts->tangent1);
1141 auto tangent_2 = getFTensor1FromMat<3>(dataAtIntegrationPts->tangent2);
1142
1143 for (unsigned int gg = 0; gg != ngp; ++gg) {
1144 normal_original_slave(i) =
1145 FTensor::levi_civita(i, j, k) * tangent_1(j) * tangent_2(k);
1146 ++normal_original_slave;
1147 ++tangent_1;
1148 ++tangent_2;
1149 }
1150
1152}
1153
1155 int row_side, int col_side, EntityType row_type, EntityType col_type,
1156 EntData &row_data, EntData &col_data) {
1157
1159
1160 if (col_type != MBVERTEX)
1162
1163 dataAtSpringPts->faceRowData = &row_data;
1164 CHKERR loopSideVolumes(sideFeName, *sideFe);
1165
1167}
1168
1171 const std::string field_name,
1172 const std::string mesh_nodals_positions) {
1174
1175 // Define boundary element that operates on rows, columns and data of a
1176 // given field
1177 CHKERR m_field.add_finite_element("SPRING", MF_ZERO);
1181 if (m_field.check_field(mesh_nodals_positions)) {
1183 mesh_nodals_positions);
1184 }
1185 // Add entities to that element, here we add all triangles with SPRING_BC
1186 // from cubit
1188 if (bit->getName().compare(0, 9, "SPRING_BC") == 0) {
1189 CHKERR m_field.add_ents_to_finite_element_by_dim(bit->getMeshset(), 2,
1190 "SPRING");
1191 }
1192 }
1193
1195}
1196
1198 MoFEM::Interface &m_field, const std::string field_name,
1199 const std::string mesh_nodals_positions, Range &spring_triangles) {
1201
1202 // Define boundary element that operates on rows, columns and data of a
1203 // given field
1204 CHKERR m_field.add_finite_element("SPRING_ALE", MF_ZERO);
1208 CHKERR m_field.modify_finite_element_add_field_row("SPRING_ALE",
1209 mesh_nodals_positions);
1210 CHKERR m_field.modify_finite_element_add_field_col("SPRING_ALE",
1211 mesh_nodals_positions);
1212 CHKERR m_field.modify_finite_element_add_field_data("SPRING_ALE",
1213 mesh_nodals_positions);
1214
1215 CHKERR m_field.add_ents_to_finite_element_by_dim(spring_triangles, 2,
1216 "SPRING_ALE");
1217
1219}
1220
1222 MoFEM::Interface &m_field,
1223 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr,
1224 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr,
1225 const std::string field_name, const std::string mesh_nodals_positions, double stiffness_scale) {
1227
1228 // Push operators to instances for springs
1229 // loop over blocks
1230 boost::shared_ptr<MetaSpringBC::DataAtIntegrationPtsSprings> commonDataPtr =
1231 boost::make_shared<MetaSpringBC::DataAtIntegrationPtsSprings>(m_field, stiffness_scale);
1232 CHKERR commonDataPtr->getParameters();
1233
1234 if (m_field.check_field(mesh_nodals_positions)) {
1235 CHKERR addHOOpsFace3D(mesh_nodals_positions, *fe_spring_lhs_ptr, false,
1236 false);
1237 CHKERR addHOOpsFace3D(mesh_nodals_positions, *fe_spring_rhs_ptr, false,
1238 false);
1239 }
1240
1241 for (auto &sitSpring : commonDataPtr->mapSpring) {
1242
1243 fe_spring_lhs_ptr->getOpPtrVector().push_back(
1244 new OpGetTangentSpEle(mesh_nodals_positions, commonDataPtr));
1245 fe_spring_lhs_ptr->getOpPtrVector().push_back(
1246 new OpGetNormalSpEle(mesh_nodals_positions, commonDataPtr));
1247 fe_spring_lhs_ptr->getOpPtrVector().push_back(
1248 new OpSpringKs(commonDataPtr, sitSpring.second, field_name));
1249
1250 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1251 new OpCalculateVectorFieldValues<3>(field_name, commonDataPtr->xAtPts));
1252 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1253 new OpCalculateVectorFieldValues<3>(mesh_nodals_positions,
1254 commonDataPtr->xInitAtPts));
1255 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1256 new OpGetTangentSpEle(mesh_nodals_positions, commonDataPtr));
1257 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1258 new OpGetNormalSpEle(mesh_nodals_positions, commonDataPtr));
1259 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1260 new OpSpringFs(commonDataPtr, sitSpring.second, field_name));
1261 }
1262
1264}
1265
1267 MoFEM::Interface &m_field,
1268 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr_dx,
1269 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_lhs_ptr_dX,
1270 boost::shared_ptr<FaceElementForcesAndSourcesCore> fe_spring_rhs_ptr,
1271 boost::shared_ptr<MetaSpringBC::DataAtIntegrationPtsSprings>
1272 data_at_integration_pts,
1273 const std::string field_name, const std::string mesh_nodals_positions,
1274 std::string side_fe_name) {
1276
1277 // Push operators to instances for springs
1278 // loop over blocks
1279 CHKERR data_at_integration_pts->getParameters();
1280
1281 if (m_field.check_field(mesh_nodals_positions)) {
1282 CHKERR addHOOpsFace3D(mesh_nodals_positions, *fe_spring_lhs_ptr_dx, false,
1283 false);
1284 CHKERR addHOOpsFace3D(mesh_nodals_positions, *fe_spring_lhs_ptr_dX, false,
1285 false);
1286 CHKERR addHOOpsFace3D(mesh_nodals_positions, *fe_spring_rhs_ptr, false,
1287 false);
1288 }
1289
1290 for (auto &sitSpring : data_at_integration_pts->mapSpring) {
1291
1292 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> feMatSideRhs =
1293 boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(m_field);
1294
1295 feMatSideRhs->getOpPtrVector().push_back(
1297 mesh_nodals_positions, data_at_integration_pts->HMat));
1298 feMatSideRhs->getOpPtrVector().push_back(
1300 field_name, data_at_integration_pts->hMat));
1301
1302 feMatSideRhs->getOpPtrVector().push_back(new OpCalculateDeformation(
1303 mesh_nodals_positions, data_at_integration_pts));
1304
1305 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1306 new OpGetTangentSpEle(mesh_nodals_positions, data_at_integration_pts));
1307 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1308 new OpGetNormalSpEle(mesh_nodals_positions, data_at_integration_pts));
1309
1310 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1312 data_at_integration_pts->xAtPts));
1313 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1315 mesh_nodals_positions, data_at_integration_pts->xInitAtPts));
1316
1317 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1318 new OpSpringFsMaterial(mesh_nodals_positions, data_at_integration_pts,
1319 feMatSideRhs, side_fe_name, sitSpring.second));
1320
1321 fe_spring_lhs_ptr_dx->getOpPtrVector().push_back(
1323 data_at_integration_pts->xAtPts));
1324 fe_spring_lhs_ptr_dx->getOpPtrVector().push_back(
1326 mesh_nodals_positions, data_at_integration_pts->xInitAtPts));
1327 fe_spring_lhs_ptr_dx->getOpPtrVector().push_back(
1328 new OpGetTangentSpEle(mesh_nodals_positions, data_at_integration_pts));
1329 fe_spring_lhs_ptr_dx->getOpPtrVector().push_back(
1330 new OpGetNormalSpEle(mesh_nodals_positions, data_at_integration_pts));
1331
1332 fe_spring_lhs_ptr_dx->getOpPtrVector().push_back(
1333 new OpSpringKs_dX(data_at_integration_pts, sitSpring.second, field_name,
1334 mesh_nodals_positions));
1335
1336 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> feMatSideLhs_dx =
1337 boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(m_field);
1338
1339 feMatSideLhs_dx->getOpPtrVector().push_back(
1341 mesh_nodals_positions, data_at_integration_pts->HMat));
1342 feMatSideLhs_dx->getOpPtrVector().push_back(
1344 field_name, data_at_integration_pts->hMat));
1345
1346 feMatSideLhs_dx->getOpPtrVector().push_back(new OpCalculateDeformation(
1347 mesh_nodals_positions, data_at_integration_pts));
1348
1349 feMatSideLhs_dx->getOpPtrVector().push_back(
1351 mesh_nodals_positions, field_name, data_at_integration_pts));
1352
1353 fe_spring_lhs_ptr_dX->getOpPtrVector().push_back(
1355 data_at_integration_pts->xAtPts));
1356 fe_spring_lhs_ptr_dX->getOpPtrVector().push_back(
1358 mesh_nodals_positions, data_at_integration_pts->xInitAtPts));
1359 fe_spring_lhs_ptr_dX->getOpPtrVector().push_back(
1360 new OpGetTangentSpEle(mesh_nodals_positions, data_at_integration_pts));
1361 fe_spring_lhs_ptr_dX->getOpPtrVector().push_back(
1362 new OpGetNormalSpEle(mesh_nodals_positions, data_at_integration_pts));
1363
1364 fe_spring_lhs_ptr_dX->getOpPtrVector().push_back(
1365 new OpSpringALEMaterialLhs_dX_dx(mesh_nodals_positions, field_name,
1366 data_at_integration_pts,
1367 feMatSideLhs_dx, side_fe_name));
1368
1369 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> feMatSideLhs_dX =
1370 boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(m_field);
1371
1372 feMatSideLhs_dX->getOpPtrVector().push_back(
1374 mesh_nodals_positions, data_at_integration_pts->HMat));
1375 feMatSideLhs_dX->getOpPtrVector().push_back(
1377 field_name, data_at_integration_pts->hMat));
1378
1379 feMatSideLhs_dX->getOpPtrVector().push_back(new OpCalculateDeformation(
1380 mesh_nodals_positions, data_at_integration_pts));
1381
1382 feMatSideLhs_dX->getOpPtrVector().push_back(
1383 new SpringALEMaterialVolOnSideLhs_dX_dX(mesh_nodals_positions,
1384 mesh_nodals_positions,
1385 data_at_integration_pts));
1386
1387 fe_spring_lhs_ptr_dX->getOpPtrVector().push_back(
1389 mesh_nodals_positions, mesh_nodals_positions,
1390 data_at_integration_pts, feMatSideLhs_dX, side_fe_name));
1391 }
1392
1394}
Kronecker Delta class symmetric.
Kronecker Delta class.
@ MF_ZERO
Definition: definitions.h:111
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ BLOCKSET
Definition: definitions.h:161
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
constexpr auto t_kd
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_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_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
#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
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
const FTensor::Tensor2< T, Dim, Dim > Vec
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
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
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:1218
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:149
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1207
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
double w(const double g, const double t)
const double r
rate factor
constexpr auto field_name
FTensor::Index< 'm', 3 > m
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