v0.13.2
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
13MoFEMErrorCode MetaSpringBC::OpSpringFs::doWork(int side, EntityType type,
14 EntData &data) {
15
16 MoFEMFunctionBegin;
17 // check that the faces have associated degrees of freedom
18 const int nb_dofs = data.getIndices().size();
19 if (nb_dofs == 0)
20 MoFEMFunctionReturnHot(0);
21 if (dAta.tRis.find(getFEEntityHandle()) == dAta.tRis.end()) {
22 MoFEMFunctionReturnHot(0);
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
36 auto t_w = getFTensor0IntegrationWeight();
37
38 // FTensor indices
39 FTensor::Index<'i', 3> i;
40 FTensor::Index<'j', 3> j;
41 FTensor::Index<'k', 3> k;
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
47 FTensor::Tensor1<double, 3> t_normal;
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
87 FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_nf(&nF[0], &nF[1],
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);
113 MoFEMFunctionReturn(0);
114}
115
116MoFEMErrorCode MetaSpringBC::OpSpringKs::doWork(int row_side, int col_side,
117 EntityType row_type,
118 EntityType col_type,
119 EntData &row_data,
120 EntData &col_data) {
121 MoFEMFunctionBegin;
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)
125 MoFEMFunctionReturnHot(0);
126
127 const int col_nb_dofs = col_data.getIndices().size();
128 if (!col_nb_dofs)
129 MoFEMFunctionReturnHot(0);
130
131 if (dAta.tRis.find(getFEEntityHandle()) == dAta.tRis.end()) {
132 MoFEMFunctionReturnHot(0);
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
143 MoFEMFunctionReturnHot(0);
144
145 // get intergration weights
146 auto t_w = getFTensor0IntegrationWeight();
147
148 FTensor::Index<'i', 3> i;
149 FTensor::Index<'j', 3> j;
150 FTensor::Index<'k', 3> k;
151
152 // create a 3d vector to be used as the normal to the face with length equal
153 // to the face area
154 FTensor::Tensor1<double, 3> t_normal;
155
156 // First tangent vector
157 auto t_tangent1_ptr = getFTensor1Tangent1AtGaussPts();
158 FTensor::Tensor1<double, 3> t_tangent1;
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
162 FTensor::Tensor1<double, 3> t_tangent2;
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
225 MoFEMFunctionReturn(0);
226}
227
228MoFEMErrorCode MetaSpringBC::OpSpringKs_dX::doWork(int row_side, int col_side,
229 EntityType row_type,
230 EntityType col_type,
231 EntData &row_data,
232 EntData &col_data) {
233 MoFEMFunctionBegin;
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)
238 MoFEMFunctionReturnHot(0);
239
240 const int col_nb_dofs = col_data.getIndices().size();
241 if (!col_nb_dofs)
242 MoFEMFunctionReturnHot(0);
243
244 if (dAta.tRis.find(getFEEntityHandle()) == dAta.tRis.end()) {
245 MoFEMFunctionReturnHot(0);
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
256 MoFEMFunctionReturnHot(0);
257
258 // get intergration weights
259 auto t_w = getFTensor0IntegrationWeight();
260
261 FTensor::Index<'i', 3> i;
262 FTensor::Index<'j', 3> j;
263 FTensor::Index<'k', 3> k;
264 FTensor::Index<'l', 3> l;
265
266 auto make_vec_der = [&](auto t_N, auto t_1, auto t_2) {
267 FTensor::Tensor2<double, 3, 3> t_n;
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) {
275 FTensor::Tensor1<double, 3> t_normal;
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));
278 FTensor::Tensor2<double, 3, 3> t_n;
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);
282 FTensor::Tensor2<double, 3, 3> t_result;
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
295 FTensor::Tensor1<double, 3> t_normal;
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
377 MoFEMFunctionReturn(0);
378}
379
380MoFEMErrorCode MetaSpringBC::SpringALEMaterialVolOnSideLhs::doWork(
381 int row_side, int col_side, EntityType row_type, EntityType col_type,
382 EntData &row_data, EntData &col_data) {
383
384 MoFEMFunctionBegin;
385
386 if (dataAtSpringPts->faceRowData == nullptr)
387 MoFEMFunctionReturnHot(0);
388
389 if (row_type != MBVERTEX)
390 MoFEMFunctionReturnHot(0);
391
392 row_nb_dofs = dataAtSpringPts->faceRowData->getIndices().size();
393 if (!row_nb_dofs)
394 MoFEMFunctionReturnHot(0);
395 col_nb_dofs = col_data.getIndices().size();
396 if (!col_nb_dofs)
397 MoFEMFunctionReturnHot(0);
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
413 MoFEMFunctionReturn(0);
414}
415
416MoFEMErrorCode
417MetaSpringBC::SpringALEMaterialVolOnSideLhs::aSsemble(EntData &row_data,
418 EntData &col_data) {
419
420 MoFEMFunctionBegin;
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
447 MoFEMFunctionReturn(0);
448}
449
450MoFEMErrorCode MetaSpringBC::SpringALEMaterialVolOnSideLhs_dX_dx::iNtegrate(
451 EntData &row_data, EntData &col_data) {
452
453 MoFEMFunctionBegin;
454
455 FTensor::Index<'i', 3> i;
456 FTensor::Index<'j', 3> j;
457 FTensor::Index<'k', 3> k;
458 FTensor::Index<'l', 3> l;
459
460 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
461 return FTensor::Tensor2<double *, 3, 3>(
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 };
466 FTensor::Tensor1<double, 3> t_normal;
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
539 MoFEMFunctionReturn(0);
540}
541
542MoFEMErrorCode
543MetaSpringBC::OpSpringALEMaterialLhs::aSsemble(EntData &row_data,
544 EntData &col_data) {
545
546 MoFEMFunctionBegin;
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
573 MoFEMFunctionReturn(0);
574}
575
576MoFEMErrorCode MetaSpringBC::OpSpringALEMaterialLhs_dX_dX::doWork(
577 int row_side, int col_side, EntityType row_type, EntityType col_type,
578 EntData &row_data, EntData &col_data) {
579
580 MoFEMFunctionBegin;
581
582 row_nb_dofs = row_data.getIndices().size();
583 if (!row_nb_dofs)
584 MoFEMFunctionReturnHot(0);
585 col_nb_dofs = col_data.getIndices().size();
586 if (!col_nb_dofs)
587 MoFEMFunctionReturnHot(0);
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
609 MoFEMFunctionReturn(0);
610}
611
612MoFEMErrorCode
613MetaSpringBC::OpSpringALEMaterialLhs_dX_dX::iNtegrate(EntData &row_data,
614 EntData &col_data) {
615
616 MoFEMFunctionBegin;
617
618 FTensor::Index<'i', 3> i;
619 FTensor::Index<'j', 3> j;
620 FTensor::Index<'k', 3> k;
621 FTensor::Index<'l', 3> l;
622
623 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
624 return FTensor::Tensor2<double *, 3, 3>(
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) {
631 FTensor::Tensor2<double, 3, 3> t_n;
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) {
639 FTensor::Tensor1<double, 3> t_normal;
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));
642 FTensor::Tensor2<double, 3, 3> t_n;
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);
646 FTensor::Tensor2<double, 3, 3> t_result;
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
668 FTensor::Tensor1<double, 3> t_normal;
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);
680 constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
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
746 MoFEMFunctionReturn(0);
747}
748
749MoFEMErrorCode MetaSpringBC::SpringALEMaterialVolOnSideLhs_dX_dX::iNtegrate(
750 EntData &row_data, EntData &col_data) {
751
752 MoFEMFunctionBegin;
753
754 FTensor::Index<'i', 3> i;
755 FTensor::Index<'j', 3> j;
756 FTensor::Index<'k', 3> k;
757 FTensor::Index<'l', 3> l;
758 FTensor::Index<'m', 3> m;
759 FTensor::Index<'q', 3> q;
760
761 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
762 return FTensor::Tensor2<double *, 3, 3>(
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
773 FTensor::Tensor1<double, 3> t_normal;
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
791 FTensor::Tensor2<double, 3, 3> t_d;
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
843 MoFEMFunctionReturn(0);
844}
845
846MoFEMErrorCode MetaSpringBC::OpCalculateDeformation::doWork(int side,
847 EntityType type,
848 EntData &row_data) {
849
850 MoFEMFunctionBegin;
851 FTensor::Index<'i', 3> i;
852 FTensor::Index<'j', 3> j;
853 FTensor::Index<'k', 3> k;
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
884 MoFEMFunctionReturn(0);
885}
886
887MoFEMErrorCode MetaSpringBC::OpSpringFsMaterial::doWork(int side,
888 EntityType type,
889 EntData &row_data) {
890
891 MoFEMFunctionBegin;
892
893 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
894 dAta.tRis.end()) {
895 MoFEMFunctionReturnHot(0);
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)
902 MoFEMFunctionReturnHot(0);
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
916 MoFEMFunctionReturn(0);
917}
918
919MoFEMErrorCode MetaSpringBC::OpSpringFsMaterial::iNtegrate(EntData &data) {
920 MoFEMFunctionBegin;
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)
927 MoFEMFunctionReturnHot(0);
928
929 CHKERR dataAtPts->getBlockData(dAta);
930
931 // get integration weights
932 auto t_w = getFTensor0IntegrationWeight();
933
934 // FTensor indices
935 FTensor::Index<'i', 3> i;
936 FTensor::Index<'j', 3> j;
937 FTensor::Index<'k', 3> k;
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
944 FTensor::Tensor1<double, 3> t_normal;
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
979 FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_nf(&nF[0], &nF[1],
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);
1006 MoFEMFunctionReturn(0);
1007}
1008
1009MoFEMErrorCode MetaSpringBC::OpSpringFsMaterial::aSsemble(EntData &row_data) {
1010 MoFEMFunctionBegin;
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
1033 MoFEMFunctionReturn(0);
1034}
1035
1036MoFEMErrorCode MetaSpringBC::OpGetTangentSpEle::doWork(int side,
1037 EntityType type,
1038 EntData &data) {
1039 MoFEMFunctionBegin;
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;
1047 FTensor::Index<'i', 3> i;
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
1073 MoFEMFunctionReturn(0);
1074}
1075
1076MoFEMErrorCode MetaSpringBC::OpGetNormalSpEle::doWork(int side, EntityType type,
1077 EntData &data) {
1078 MoFEMFunctionBegin;
1079
1080 if (data.getFieldData().size() == 0)
1081 MoFEMFunctionReturnHot(0);
1082
1083 ngp = data.getN().size1();
1084
1085 FTensor::Index<'i', 3> i;
1086 FTensor::Index<'j', 3> j;
1087 FTensor::Index<'k', 3> k;
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
1108 MoFEMFunctionReturn(0);
1109}
1110
1111MoFEMErrorCode MetaSpringBC::OpSpringALEMaterialLhs_dX_dx::doWork(
1112 int row_side, int col_side, EntityType row_type, EntityType col_type,
1113 EntData &row_data, EntData &col_data) {
1114
1115 MoFEMFunctionBegin;
1116
1117 if (col_type != MBVERTEX)
1118 MoFEMFunctionReturnHot(0);
1119
1120 dataAtSpringPts->faceRowData = &row_data;
1121 CHKERR loopSideVolumes(sideFeName, *sideFe);
1122
1123 MoFEMFunctionReturn(0);
1124}
1125
1126MoFEMErrorCode
1127MetaSpringBC::addSpringElements(MoFEM::Interface &m_field,
1128 const std::string field_name,
1129 const std::string mesh_nodals_positions) {
1130 MoFEMFunctionBegin;
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);
1135 CHKERR m_field.modify_finite_element_add_field_row("SPRING", field_name);
1136 CHKERR m_field.modify_finite_element_add_field_col("SPRING", field_name);
1137 CHKERR m_field.modify_finite_element_add_field_data("SPRING", field_name);
1138 if (m_field.check_field(mesh_nodals_positions)) {
1139 CHKERR m_field.modify_finite_element_add_field_data("SPRING",
1140 mesh_nodals_positions);
1141 }
1142 // Add entities to that element, here we add all triangles with SPRING_BC
1143 // from cubit
1144 for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, bit)) {
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
1151 MoFEMFunctionReturn(0);
1152}
1153
1154MoFEMErrorCode MetaSpringBC::addSpringElementsALE(
1155 MoFEM::Interface &m_field, const std::string field_name,
1156 const std::string mesh_nodals_positions, Range &spring_triangles) {
1157 MoFEMFunctionBegin;
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);
1162 CHKERR m_field.modify_finite_element_add_field_row("SPRING_ALE", field_name);
1163 CHKERR m_field.modify_finite_element_add_field_col("SPRING_ALE", field_name);
1164 CHKERR m_field.modify_finite_element_add_field_data("SPRING_ALE", field_name);
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
1175 MoFEMFunctionReturn(0);
1176}
1177
1178MoFEMErrorCode MetaSpringBC::setSpringOperators(
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) {
1183 MoFEMFunctionBegin;
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
1220 MoFEMFunctionReturn(0);
1221}
1222
1223MoFEMErrorCode MetaSpringBC::setSpringOperatorsMaterial(
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) {
1232 MoFEMFunctionBegin;
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(
1253 new OpCalculateVectorFieldGradient<3, 3>(
1254 mesh_nodals_positions, data_at_integration_pts->HMat));
1255 feMatSideRhs->getOpPtrVector().push_back(
1256 new OpCalculateVectorFieldGradient<3, 3>(
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(
1268 new OpCalculateVectorFieldValues<3>(field_name,
1269 data_at_integration_pts->xAtPts));
1270 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1271 new OpCalculateVectorFieldValues<3>(
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(
1279 new OpCalculateVectorFieldValues<3>(field_name,
1280 data_at_integration_pts->xAtPts));
1281 fe_spring_lhs_ptr_dx->getOpPtrVector().push_back(
1282 new OpCalculateVectorFieldValues<3>(
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(
1297 new OpCalculateVectorFieldGradient<3, 3>(
1298 mesh_nodals_positions, data_at_integration_pts->HMat));
1299 feMatSideLhs_dx->getOpPtrVector().push_back(
1300 new OpCalculateVectorFieldGradient<3, 3>(
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(
1307 new SpringALEMaterialVolOnSideLhs_dX_dx(
1308 mesh_nodals_positions, field_name, data_at_integration_pts));
1309
1310 fe_spring_lhs_ptr_dX->getOpPtrVector().push_back(
1311 new OpCalculateVectorFieldValues<3>(field_name,
1312 data_at_integration_pts->xAtPts));
1313 fe_spring_lhs_ptr_dX->getOpPtrVector().push_back(
1314 new OpCalculateVectorFieldValues<3>(
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(
1330 new OpCalculateVectorFieldGradient<3, 3>(
1331 mesh_nodals_positions, data_at_integration_pts->HMat));
1332 feMatSideLhs_dX->getOpPtrVector().push_back(
1333 new OpCalculateVectorFieldGradient<3, 3>(
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(
1345 new OpSpringALEMaterialLhs_dX_dX(
1346 mesh_nodals_positions, mesh_nodals_positions,
1347 data_at_integration_pts, feMatSideLhs_dX, side_fe_name));
1348 }
1349
1350 MoFEMFunctionReturn(0);
1351}
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
FTensor::Kronecker_Delta
Kronecker Delta class.
Definition: Kronecker_Delta.hpp:15
MF_ZERO
@ MF_ZERO
Definition: definitions.h:98
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
BLOCKSET
@ BLOCKSET
Definition: definitions.h:148
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
m
FTensor::Index< 'm', SPACE_DIM > m
Definition: fluid_structure_eigenproblem.cpp:65
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:73
MoFEM::CoreInterface::add_ents_to_finite_element_by_dim
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
MoFEM::CoreInterface::add_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
MoFEM::CoreInterface::modify_finite_element_add_field_col
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
MoFEM::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
MoFEM::CoreInterface::modify_finite_element_add_field_row
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
MoFEM::CoreInterface::check_field
virtual bool check_field(const std::string &name) const =0
check if field is in database
_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_
#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.
Definition: MeshsetsManager.hpp:71
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
FTensor::levi_civita
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
Definition: Levi_Civita.hpp:617
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::addHOOpsFace3D
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
Definition: HODataOperators.hpp:754
MoFEM::invertTensor3by3
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:1363
MoFEM::MatSetValues
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
Definition: EntitiesFieldData.hpp:1631
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1352
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1576
MoFEM::VectorDofs
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
Definition: EntitiesFieldData.hpp:23
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
MetaSpringBC::OpCalculateDeformation
Operator for computing deformation gradients in side volumes.
Definition: SpringElement.hpp:681
MetaSpringBC::OpCalculateDeformation::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &row_data)
Definition: SpringElement.cpp:846
MetaSpringBC::OpGetNormalSpEle
Computes, for material configuration, normal to face that lies on a surface with springs.
Definition: SpringElement.hpp:861
MetaSpringBC::OpGetNormalSpEle::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &row_data)
Evaluates normal vector to the triangle on surface with springs based on material base coordinates.
Definition: SpringElement.cpp:1076
MetaSpringBC::OpGetTangentSpEle
Computes, for material configuration, tangent vectors to face that lie on a surface with springs.
Definition: SpringElement.hpp:802
MetaSpringBC::OpGetTangentSpEle::doWork
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...
Definition: SpringElement.cpp:1036
MetaSpringBC::OpSpringALEMaterialLhs_dX_dX
LHS-operator for the pressure element (material configuration)
Definition: SpringElement.hpp:544
MetaSpringBC::OpSpringALEMaterialLhs_dX_dX::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Definition: SpringElement.cpp:576
MetaSpringBC::OpSpringALEMaterialLhs_dX_dX::iNtegrate
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Compute part of the left-hand side.
Definition: SpringElement.cpp:613
MetaSpringBC::OpSpringALEMaterialLhs_dX_dx
LHS-operator for the pressure element (material configuration)
Definition: SpringElement.hpp:515
MetaSpringBC::OpSpringALEMaterialLhs_dX_dx::doWork
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.
Definition: SpringElement.cpp:1111
MetaSpringBC::OpSpringALEMaterialLhs::aSsemble
MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data)
Definition: SpringElement.cpp:543
MetaSpringBC::OpSpringFs
RHS-operator for the spring boundary condition element.
Definition: SpringElement.hpp:134
MetaSpringBC::OpSpringFs::dAta
MetaSpringBC::BlockOptionDataSprings & dAta
Definition: SpringElement.hpp:140
MetaSpringBC::OpSpringFs::commonDataPtr
boost::shared_ptr< MetaSpringBC::DataAtIntegrationPtsSprings > commonDataPtr
Definition: SpringElement.hpp:139
MetaSpringBC::OpSpringFs::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Integrates and assembles to global RHS vector virtual work of springs.
Definition: SpringElement.cpp:13
MetaSpringBC::OpSpringFsMaterial
RHS-operator for the spring boundary condition element for ALE formulation.
Definition: SpringElement.hpp:716
MetaSpringBC::OpSpringFsMaterial::aSsemble
MoFEMErrorCode aSsemble(EntData &row_data)
Definition: SpringElement.cpp:1009
MetaSpringBC::OpSpringFsMaterial::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
Definition: SpringElement.cpp:887
MetaSpringBC::OpSpringFsMaterial::iNtegrate
MoFEMErrorCode iNtegrate(EntData &row_data)
Integrates and assembles to global RHS vector virtual work of springs on material positions for confi...
Definition: SpringElement.cpp:919
MetaSpringBC::OpSpringKs_dX
Definition: SpringElement.hpp:255
MetaSpringBC::OpSpringKs_dX::doWork
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.
Definition: SpringElement.cpp:228
MetaSpringBC::OpSpringKs
LHS-operator for the springs element.
Definition: SpringElement.hpp:200
MetaSpringBC::OpSpringKs::doWork
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.
Definition: SpringElement.cpp:116
MetaSpringBC::SpringALEMaterialVolOnSideLhs_dX_dX
LHS-operator (material configuration) on the side volume.
Definition: SpringElement.hpp:631
MetaSpringBC::SpringALEMaterialVolOnSideLhs_dX_dX::iNtegrate
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrates over a face contribution from a side volume.
Definition: SpringElement.cpp:749
MetaSpringBC::SpringALEMaterialVolOnSideLhs_dX_dx
LHS-operator (material configuration) on the side volume for spring element.
Definition: SpringElement.hpp:389
MetaSpringBC::SpringALEMaterialVolOnSideLhs_dX_dx::iNtegrate
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrates over a face contribution from a side volume.
Definition: SpringElement.cpp:450
MetaSpringBC::SpringALEMaterialVolOnSideLhs::aSsemble
MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data)
Definition: SpringElement.cpp:417
MetaSpringBC::SpringALEMaterialVolOnSideLhs::doWork
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.
Definition: SpringElement.cpp:380
MetaSpringBC::setSpringOperators
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.
Definition: SpringElement.cpp:1178
MetaSpringBC::addSpringElementsALE
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.
Definition: SpringElement.cpp:1154
MetaSpringBC::addSpringElements
static MoFEMErrorCode addSpringElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Declare spring element.
Definition: SpringElement.cpp:1127
MetaSpringBC::setSpringOperatorsMaterial
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.
Definition: SpringElement.cpp:1223
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:127
MoFEM::EntitiesFieldData::EntData::getFTensor1DiffN
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
Definition: EntitiesFieldData.cpp:387
MoFEM::EntitiesFieldData::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition: EntitiesFieldData.hpp:1489
MoFEM::EntitiesFieldData::EntData::getN
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
Definition: EntitiesFieldData.hpp:1305
MoFEM::EntitiesFieldData::EntData::getFieldData
const VectorDouble & getFieldData() const
get dofs values
Definition: EntitiesFieldData.hpp:1241
MoFEM::EntitiesFieldData::EntData::getFTensor1FieldData
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1FieldData()
Return FTensor of rank 1, i.e. vector from filed data coeffinects.
Definition: EntitiesFieldData.hpp:1275
MoFEM::EntitiesFieldData::EntData::getFieldDofs
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
Definition: EntitiesFieldData.hpp:1256
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition: EntitiesFieldData.hpp:1201
MoFEM::ForcesAndSourcesCore::UserDataOperator::getFEEntityHandle
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
Definition: ForcesAndSourcesCore.hpp:992
MoFEM::ForcesAndSourcesCore::UserDataOperator::getFTensor0IntegrationWeight
auto getFTensor0IntegrationWeight()
Get integration weights.
Definition: ForcesAndSourcesCore.hpp:1228
MoFEM::ForcesAndSourcesCore::UserDataOperator::getNumeredEntFiniteElementPtr
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
Definition: ForcesAndSourcesCore.hpp:988
MoFEM::ForcesAndSourcesCore::UserDataOperator::getFEMethod
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
Definition: ForcesAndSourcesCore.hpp:1030
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1513
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:450
MoFEM::SnesMethod::snes_f
Vec & snes_f
residual
Definition: LoopMethods.hpp:122
Generated by Doxygen 1.9.5 and hosted at University of Glasgow