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
25MoFEMErrorCode MetaSpringBC::OpSpringFs::doWork(int side, EntityType type,
26 EntData &data) {
27
28 MoFEMFunctionBegin;
29 // check that the faces have associated degrees of freedom
30 const int nb_dofs = data.getIndices().size();
31 if (nb_dofs == 0)
32 MoFEMFunctionReturnHot(0);
33 if (dAta.tRis.find(getFEEntityHandle()) == dAta.tRis.end()) {
34 MoFEMFunctionReturnHot(0);
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
48 auto t_w = getFTensor0IntegrationWeight();
49
50 // FTensor indices
51 FTensor::Index<'i', 3> i;
52 FTensor::Index<'j', 3> j;
53 FTensor::Index<'k', 3> k;
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
59 FTensor::Tensor1<double, 3> t_normal;
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
99 FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_nf(&nF[0], &nF[1],
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);
125 MoFEMFunctionReturn(0);
126}
127
128MoFEMErrorCode MetaSpringBC::OpSpringKs::doWork(int row_side, int col_side,
129 EntityType row_type,
130 EntityType col_type,
131 EntData &row_data,
132 EntData &col_data) {
133 MoFEMFunctionBegin;
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)
137 MoFEMFunctionReturnHot(0);
138
139 const int col_nb_dofs = col_data.getIndices().size();
140 if (!col_nb_dofs)
141 MoFEMFunctionReturnHot(0);
142
143 if (dAta.tRis.find(getFEEntityHandle()) == dAta.tRis.end()) {
144 MoFEMFunctionReturnHot(0);
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
155 MoFEMFunctionReturnHot(0);
156
157 // get intergration weights
158 auto t_w = getFTensor0IntegrationWeight();
159
160 FTensor::Index<'i', 3> i;
161 FTensor::Index<'j', 3> j;
162 FTensor::Index<'k', 3> k;
163
164 // create a 3d vector to be used as the normal to the face with length equal
165 // to the face area
166 FTensor::Tensor1<double, 3> t_normal;
167
168 // First tangent vector
169 auto t_tangent1_ptr = getFTensor1Tangent1AtGaussPts();
170 FTensor::Tensor1<double, 3> t_tangent1;
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
174 FTensor::Tensor1<double, 3> t_tangent2;
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
237 MoFEMFunctionReturn(0);
238}
239
240MoFEMErrorCode MetaSpringBC::OpSpringKs_dX::doWork(int row_side, int col_side,
241 EntityType row_type,
242 EntityType col_type,
243 EntData &row_data,
244 EntData &col_data) {
245 MoFEMFunctionBegin;
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)
250 MoFEMFunctionReturnHot(0);
251
252 const int col_nb_dofs = col_data.getIndices().size();
253 if (!col_nb_dofs)
254 MoFEMFunctionReturnHot(0);
255
256 if (dAta.tRis.find(getFEEntityHandle()) == dAta.tRis.end()) {
257 MoFEMFunctionReturnHot(0);
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
268 MoFEMFunctionReturnHot(0);
269
270 // get intergration weights
271 auto t_w = getFTensor0IntegrationWeight();
272
273 FTensor::Index<'i', 3> i;
274 FTensor::Index<'j', 3> j;
275 FTensor::Index<'k', 3> k;
276 FTensor::Index<'l', 3> l;
277
278 auto make_vec_der = [&](auto t_N, auto t_1, auto t_2) {
279 FTensor::Tensor2<double, 3, 3> t_n;
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) {
287 FTensor::Tensor1<double, 3> t_normal;
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));
290 FTensor::Tensor2<double, 3, 3> t_n;
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);
294 FTensor::Tensor2<double, 3, 3> t_result;
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
307 FTensor::Tensor1<double, 3> t_normal;
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
389 MoFEMFunctionReturn(0);
390}
391
392MoFEMErrorCode MetaSpringBC::SpringALEMaterialVolOnSideLhs::doWork(
393 int row_side, int col_side, EntityType row_type, EntityType col_type,
394 EntData &row_data, EntData &col_data) {
395
396 MoFEMFunctionBegin;
397
398 if (dataAtSpringPts->faceRowData == nullptr)
399 MoFEMFunctionReturnHot(0);
400
401 if (row_type != MBVERTEX)
402 MoFEMFunctionReturnHot(0);
403
404 row_nb_dofs = dataAtSpringPts->faceRowData->getIndices().size();
405 if (!row_nb_dofs)
406 MoFEMFunctionReturnHot(0);
407 col_nb_dofs = col_data.getIndices().size();
408 if (!col_nb_dofs)
409 MoFEMFunctionReturnHot(0);
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
425 MoFEMFunctionReturn(0);
426}
427
428MoFEMErrorCode
429MetaSpringBC::SpringALEMaterialVolOnSideLhs::aSsemble(EntData &row_data,
430 EntData &col_data) {
431
432 MoFEMFunctionBegin;
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
474 MoFEMFunctionReturn(0);
475}
476
477MoFEMErrorCode MetaSpringBC::SpringALEMaterialVolOnSideLhs_dX_dx::iNtegrate(
478 EntData &row_data, EntData &col_data) {
479
480 MoFEMFunctionBegin;
481
482 FTensor::Index<'i', 3> i;
483 FTensor::Index<'j', 3> j;
484 FTensor::Index<'k', 3> k;
485 FTensor::Index<'l', 3> l;
486
487 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
488 return FTensor::Tensor2<double *, 3, 3>(
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 };
493 FTensor::Tensor1<double, 3> t_normal;
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
566 MoFEMFunctionReturn(0);
567}
568
569MoFEMErrorCode
570MetaSpringBC::OpSpringALEMaterialLhs::aSsemble(EntData &row_data,
571 EntData &col_data) {
572
573 MoFEMFunctionBegin;
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
615 MoFEMFunctionReturn(0);
616}
617
618MoFEMErrorCode MetaSpringBC::OpSpringALEMaterialLhs_dX_dX::doWork(
619 int row_side, int col_side, EntityType row_type, EntityType col_type,
620 EntData &row_data, EntData &col_data) {
621
622 MoFEMFunctionBegin;
623
624 row_nb_dofs = row_data.getIndices().size();
625 if (!row_nb_dofs)
626 MoFEMFunctionReturnHot(0);
627 col_nb_dofs = col_data.getIndices().size();
628 if (!col_nb_dofs)
629 MoFEMFunctionReturnHot(0);
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
651 MoFEMFunctionReturn(0);
652}
653
654MoFEMErrorCode
655MetaSpringBC::OpSpringALEMaterialLhs_dX_dX::iNtegrate(EntData &row_data,
656 EntData &col_data) {
657
658 MoFEMFunctionBegin;
659
660 FTensor::Index<'i', 3> i;
661 FTensor::Index<'j', 3> j;
662 FTensor::Index<'k', 3> k;
663 FTensor::Index<'l', 3> l;
664
665 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
666 return FTensor::Tensor2<double *, 3, 3>(
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) {
673 FTensor::Tensor2<double, 3, 3> t_n;
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) {
681 FTensor::Tensor1<double, 3> t_normal;
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));
684 FTensor::Tensor2<double, 3, 3> t_n;
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);
688 FTensor::Tensor2<double, 3, 3> t_result;
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
710 FTensor::Tensor1<double, 3> t_normal;
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);
722 constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
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
788 MoFEMFunctionReturn(0);
789}
790
791MoFEMErrorCode MetaSpringBC::SpringALEMaterialVolOnSideLhs_dX_dX::iNtegrate(
792 EntData &row_data, EntData &col_data) {
793
794 MoFEMFunctionBegin;
795
796 FTensor::Index<'i', 3> i;
797 FTensor::Index<'j', 3> j;
798 FTensor::Index<'k', 3> k;
799 FTensor::Index<'l', 3> l;
800 FTensor::Index<'m', 3> m;
801 FTensor::Index<'q', 3> q;
802
803 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
804 return FTensor::Tensor2<double *, 3, 3>(
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
815 FTensor::Tensor1<double, 3> t_normal;
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
833 FTensor::Tensor2<double, 3, 3> t_d;
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
885 MoFEMFunctionReturn(0);
886}
887
888MoFEMErrorCode MetaSpringBC::OpCalculateDeformation::doWork(int side,
889 EntityType type,
890 EntData &row_data) {
891
892 MoFEMFunctionBegin;
893 FTensor::Index<'i', 3> i;
894 FTensor::Index<'j', 3> j;
895 FTensor::Index<'k', 3> k;
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
926 MoFEMFunctionReturn(0);
927}
928
929MoFEMErrorCode MetaSpringBC::OpSpringFsMaterial::doWork(int side,
930 EntityType type,
931 EntData &row_data) {
932
933 MoFEMFunctionBegin;
934
935 if (dAta.tRis.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
936 dAta.tRis.end()) {
937 MoFEMFunctionReturnHot(0);
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)
944 MoFEMFunctionReturnHot(0);
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
958 MoFEMFunctionReturn(0);
959}
960
961MoFEMErrorCode MetaSpringBC::OpSpringFsMaterial::iNtegrate(EntData &data) {
962 MoFEMFunctionBegin;
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)
969 MoFEMFunctionReturnHot(0);
970
971 CHKERR dataAtPts->getBlockData(dAta);
972
973 // get integration weights
974 auto t_w = getFTensor0IntegrationWeight();
975
976 // FTensor indices
977 FTensor::Index<'i', 3> i;
978 FTensor::Index<'j', 3> j;
979 FTensor::Index<'k', 3> k;
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
986 FTensor::Tensor1<double, 3> t_normal;
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);
1048 MoFEMFunctionReturn(0);
1049}
1050
1051MoFEMErrorCode MetaSpringBC::OpSpringFsMaterial::aSsemble(EntData &row_data) {
1052 MoFEMFunctionBegin;
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
1076 MoFEMFunctionReturn(0);
1077}
1078
1079MoFEMErrorCode MetaSpringBC::OpGetTangentSpEle::doWork(int side,
1080 EntityType type,
1081 EntData &data) {
1082 MoFEMFunctionBegin;
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;
1090 FTensor::Index<'i', 3> i;
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
1116 MoFEMFunctionReturn(0);
1117}
1118
1119MoFEMErrorCode MetaSpringBC::OpGetNormalSpEle::doWork(int side, EntityType type,
1120 EntData &data) {
1121 MoFEMFunctionBegin;
1122
1123 if (data.getFieldData().size() == 0)
1124 MoFEMFunctionReturnHot(0);
1125
1126 ngp = data.getN().size1();
1127
1128 FTensor::Index<'i', 3> i;
1129 FTensor::Index<'j', 3> j;
1130 FTensor::Index<'k', 3> k;
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
1151 MoFEMFunctionReturn(0);
1152}
1153
1154MoFEMErrorCode MetaSpringBC::OpSpringALEMaterialLhs_dX_dx::doWork(
1155 int row_side, int col_side, EntityType row_type, EntityType col_type,
1156 EntData &row_data, EntData &col_data) {
1157
1158 MoFEMFunctionBegin;
1159
1160 if (col_type != MBVERTEX)
1161 MoFEMFunctionReturnHot(0);
1162
1163 dataAtSpringPts->faceRowData = &row_data;
1164 CHKERR loopSideVolumes(sideFeName, *sideFe);
1165
1166 MoFEMFunctionReturn(0);
1167}
1168
1169MoFEMErrorCode
1170MetaSpringBC::addSpringElements(MoFEM::Interface &m_field,
1171 const std::string field_name,
1172 const std::string mesh_nodals_positions) {
1173 MoFEMFunctionBegin;
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);
1178 CHKERR m_field.modify_finite_element_add_field_row("SPRING", field_name);
1179 CHKERR m_field.modify_finite_element_add_field_col("SPRING", field_name);
1180 CHKERR m_field.modify_finite_element_add_field_data("SPRING", field_name);
1181 if (m_field.check_field(mesh_nodals_positions)) {
1182 CHKERR m_field.modify_finite_element_add_field_data("SPRING",
1183 mesh_nodals_positions);
1184 }
1185 // Add entities to that element, here we add all triangles with SPRING_BC
1186 // from cubit
1187 for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, bit)) {
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
1194 MoFEMFunctionReturn(0);
1195}
1196
1197MoFEMErrorCode MetaSpringBC::addSpringElementsALE(
1198 MoFEM::Interface &m_field, const std::string field_name,
1199 const std::string mesh_nodals_positions, Range &spring_triangles) {
1200 MoFEMFunctionBegin;
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);
1205 CHKERR m_field.modify_finite_element_add_field_row("SPRING_ALE", field_name);
1206 CHKERR m_field.modify_finite_element_add_field_col("SPRING_ALE", field_name);
1207 CHKERR m_field.modify_finite_element_add_field_data("SPRING_ALE", field_name);
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
1218 MoFEMFunctionReturn(0);
1219}
1220
1221MoFEMErrorCode MetaSpringBC::setSpringOperators(
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) {
1226 MoFEMFunctionBegin;
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
1263 MoFEMFunctionReturn(0);
1264}
1265
1266MoFEMErrorCode MetaSpringBC::setSpringOperatorsMaterial(
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) {
1275 MoFEMFunctionBegin;
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(
1296 new OpCalculateVectorFieldGradient<3, 3>(
1297 mesh_nodals_positions, data_at_integration_pts->HMat));
1298 feMatSideRhs->getOpPtrVector().push_back(
1299 new OpCalculateVectorFieldGradient<3, 3>(
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(
1311 new OpCalculateVectorFieldValues<3>(field_name,
1312 data_at_integration_pts->xAtPts));
1313 fe_spring_rhs_ptr->getOpPtrVector().push_back(
1314 new OpCalculateVectorFieldValues<3>(
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(
1322 new OpCalculateVectorFieldValues<3>(field_name,
1323 data_at_integration_pts->xAtPts));
1324 fe_spring_lhs_ptr_dx->getOpPtrVector().push_back(
1325 new OpCalculateVectorFieldValues<3>(
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(
1340 new OpCalculateVectorFieldGradient<3, 3>(
1341 mesh_nodals_positions, data_at_integration_pts->HMat));
1342 feMatSideLhs_dx->getOpPtrVector().push_back(
1343 new OpCalculateVectorFieldGradient<3, 3>(
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(
1350 new SpringALEMaterialVolOnSideLhs_dX_dx(
1351 mesh_nodals_positions, field_name, data_at_integration_pts));
1352
1353 fe_spring_lhs_ptr_dX->getOpPtrVector().push_back(
1354 new OpCalculateVectorFieldValues<3>(field_name,
1355 data_at_integration_pts->xAtPts));
1356 fe_spring_lhs_ptr_dX->getOpPtrVector().push_back(
1357 new OpCalculateVectorFieldValues<3>(
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(
1373 new OpCalculateVectorFieldGradient<3, 3>(
1374 mesh_nodals_positions, data_at_integration_pts->HMat));
1375 feMatSideLhs_dX->getOpPtrVector().push_back(
1376 new OpCalculateVectorFieldGradient<3, 3>(
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(
1388 new OpSpringALEMaterialLhs_dX_dX(
1389 mesh_nodals_positions, mesh_nodals_positions,
1390 data_at_integration_pts, feMatSideLhs_dX, side_fe_name));
1391 }
1392
1393 MoFEMFunctionReturn(0);
1394}
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:111
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
BLOCKSET
@ BLOCKSET
Definition: definitions.h:161
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:90
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::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_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_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::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:81
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:102
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.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
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:78
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
FTensor::dd
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
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
MoFEM::addHOOpsFace3D
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
Definition: HODataOperators.hpp:517
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:1218
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:149
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1207
MoFEM::VectorDofs
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
Definition: EntitiesFieldData.hpp:35
OpContactTools::w
double w(const double g, const double t)
Definition: ContactOperators.hpp:109
ReactionDiffusionEquation::r
const double r
rate factor
Definition: reaction_diffusion.cpp:33
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:31
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:95
MetaSpringBC::OpCalculateDeformation
Operator for computing deformation gradients in side volumes.
Definition: SpringElement.hpp:694
MetaSpringBC::OpCalculateDeformation::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &row_data)
Definition: SpringElement.cpp:888
MetaSpringBC::OpGetNormalSpEle
Computes, for material configuration, normal to face that lies on a surface with springs.
Definition: SpringElement.hpp:874
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:1119
MetaSpringBC::OpGetTangentSpEle
Computes, for material configuration, tangent vectors to face that lie on a surface with springs.
Definition: SpringElement.hpp:815
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:1079
MetaSpringBC::OpSpringALEMaterialLhs_dX_dX
LHS-operator for the pressure element (material configuration)
Definition: SpringElement.hpp:557
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:618
MetaSpringBC::OpSpringALEMaterialLhs_dX_dX::iNtegrate
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Compute part of the left-hand side.
Definition: SpringElement.cpp:655
MetaSpringBC::OpSpringALEMaterialLhs_dX_dx
LHS-operator for the pressure element (material configuration)
Definition: SpringElement.hpp:528
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:1154
MetaSpringBC::OpSpringALEMaterialLhs::aSsemble
MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data)
Definition: SpringElement.cpp:570
MetaSpringBC::OpSpringFs
RHS-operator for the spring boundary condition element.
Definition: SpringElement.hpp:147
MetaSpringBC::OpSpringFs::dAta
MetaSpringBC::BlockOptionDataSprings & dAta
Definition: SpringElement.hpp:153
MetaSpringBC::OpSpringFs::commonDataPtr
boost::shared_ptr< MetaSpringBC::DataAtIntegrationPtsSprings > commonDataPtr
Definition: SpringElement.hpp:152
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:25
MetaSpringBC::OpSpringFsMaterial
RHS-operator for the spring boundary condition element for ALE formulation.
Definition: SpringElement.hpp:729
MetaSpringBC::OpSpringFsMaterial::aSsemble
MoFEMErrorCode aSsemble(EntData &row_data)
Definition: SpringElement.cpp:1051
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:929
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:961
MetaSpringBC::OpSpringKs_dX
Definition: SpringElement.hpp:268
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:240
MetaSpringBC::OpSpringKs
LHS-operator for the springs element.
Definition: SpringElement.hpp:213
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:128
MetaSpringBC::SpringALEMaterialVolOnSideLhs_dX_dX
LHS-operator (material configuration) on the side volume.
Definition: SpringElement.hpp:644
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:791
MetaSpringBC::SpringALEMaterialVolOnSideLhs_dX_dx
LHS-operator (material configuration) on the side volume for spring element.
Definition: SpringElement.hpp:402
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:477
MetaSpringBC::SpringALEMaterialVolOnSideLhs::aSsemble
MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data)
Definition: SpringElement.cpp:429
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:392
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:1221
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:1197
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:1170
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:1266
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:26
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:139
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:399
MoFEM::EntitiesFieldData::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition: EntitiesFieldData.hpp:1501
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:1317
MoFEM::EntitiesFieldData::EntData::getFieldData
const VectorDouble & getFieldData() const
get dofs values
Definition: EntitiesFieldData.hpp:1253
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:1287
MoFEM::EntitiesFieldData::EntData::getFieldDofs
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
Definition: EntitiesFieldData.hpp:1268
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition: EntitiesFieldData.hpp:1213
MoFEM::ForcesAndSourcesCore::UserDataOperator::getFEEntityHandle
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
Definition: ForcesAndSourcesCore.hpp:1025
MoFEM::ForcesAndSourcesCore::UserDataOperator::getFTensor0IntegrationWeight
auto getFTensor0IntegrationWeight()
Get integration weights.
Definition: ForcesAndSourcesCore.hpp:1253
MoFEM::ForcesAndSourcesCore::UserDataOperator::getNumeredEntFiniteElementPtr
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
Definition: ForcesAndSourcesCore.hpp:1021
MoFEM::ForcesAndSourcesCore::UserDataOperator::getFEMethod
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
Definition: ForcesAndSourcesCore.hpp:1063
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1460
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:422
MoFEM::SnesMethod::snes_f
Vec & snes_f
residual
Definition: LoopMethods.hpp:132
Generated on Tue Jul 5 2022 22:52:10 for MoFEM by Doxygen 1.9.4 and hosted at University of Glasgow