v0.14.0
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>
8 using namespace MoFEM;
9 
10 #include <SpringElement.hpp>
11 using namespace boost::numeric;
12 
13 MoFEMErrorCode 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 
116 MoFEMErrorCode 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 
228 MoFEMErrorCode 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 
380 MoFEMErrorCode 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 
416 MoFEMErrorCode
417 MetaSpringBC::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 
450 MoFEMErrorCode 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 
542 MoFEMErrorCode
543 MetaSpringBC::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 
576 MoFEMErrorCode 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 
612 MoFEMErrorCode
613 MetaSpringBC::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 
749 MoFEMErrorCode 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 
846 MoFEMErrorCode 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 
887 MoFEMErrorCode 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 
919 MoFEMErrorCode 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 
1009 MoFEMErrorCode 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 
1036 MoFEMErrorCode 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 
1076 MoFEMErrorCode 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 
1111 MoFEMErrorCode 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 
1126 MoFEMErrorCode
1127 MetaSpringBC::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 
1154 MoFEMErrorCode 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 
1178 MoFEMErrorCode 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 
1223 MoFEMErrorCode 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 }
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::addHOOpsFace3D
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
Definition: HODataOperators.hpp:699
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
MetaSpringBC::OpSpringFs
RHS-operator for the spring boundary condition element.
Definition: SpringElement.hpp:134
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:1644
FTensor::Tensor1< double, 3 >
MetaSpringBC::OpSpringALEMaterialLhs_dX_dX::iNtegrate
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Compute part of the left-hand side.
Definition: SpringElement.cpp:613
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
MetaSpringBC::OpSpringKs_dX
Definition: SpringElement.hpp:254
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::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
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:466
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MoFEM::EntitiesFieldData::EntData::getFieldData
const VectorDouble & getFieldData() const
get dofs values
Definition: EntitiesFieldData.hpp:1254
FTensor::Kronecker_Delta
Kronecker Delta class.
Definition: Kronecker_Delta.hpp:15
MetaSpringBC::OpSpringFsMaterial
RHS-operator for the spring boundary condition element for ALE formulation.
Definition: SpringElement.hpp:715
MoFEM.hpp
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::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
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::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1589
MetaSpringBC::OpSpringALEMaterialLhs::aSsemble
MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data)
Definition: SpringElement.cpp:543
sdf.r
int r
Definition: sdf.py:8
FTensor::Tensor2< double, 3, 3 >
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
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:1588
MoFEM::EntitiesFieldData::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition: EntitiesFieldData.hpp:1502
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MetaSpringBC::SpringALEMaterialVolOnSideLhs_dX_dX
LHS-operator (material configuration) on the side volume.
Definition: SpringElement.hpp:630
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::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
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
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
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
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
MetaSpringBC::OpGetNormalSpEle
Computes, for material configuration, normal to face that lies on a surface with springs.
Definition: SpringElement.hpp:860
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
convert.type
type
Definition: convert.py:64
MetaSpringBC::SpringALEMaterialVolOnSideLhs_dX_dx
LHS-operator (material configuration) on the side volume for spring element.
Definition: SpringElement.hpp:388
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
MetaSpringBC::SpringALEMaterialVolOnSideLhs::aSsemble
MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data)
Definition: SpringElement.cpp:417
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition: EntitiesFieldData.hpp:1214
MoFEM::VectorDofs
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
Definition: EntitiesFieldData.hpp:23
MetaSpringBC::OpSpringFsMaterial::aSsemble
MoFEMErrorCode aSsemble(EntData &row_data)
Definition: SpringElement.cpp:1009
SpringElement.hpp
MoFEM::CoreInterface::check_field
virtual bool check_field(const std::string &name) const =0
check if field is in database
MoFEM::EntitiesFieldData::EntData::getFieldDofs
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
Definition: EntitiesFieldData.hpp:1269
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::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:137
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
FTensor::Index< 'i', 3 >
MetaSpringBC::OpCalculateDeformation::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &row_data)
Definition: SpringElement.cpp:846
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1535
MetaSpringBC::OpCalculateDeformation
Operator for computing deformation gradients in side volumes.
Definition: SpringElement.hpp:679
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::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1540
Range
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
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
MF_ZERO
@ MF_ZERO
Definition: definitions.h:111
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:526
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 coefficients.
Definition: EntitiesFieldData.hpp:1288
FTensor::Tensor0
Definition: Tensor0.hpp:16
_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
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:1318
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
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
MetaSpringBC::OpGetTangentSpEle
Computes, for material configuration, tangent vectors to face that lie on a surface with springs.
Definition: SpringElement.hpp:801
BLOCKSET
@ BLOCKSET
Definition: definitions.h:161
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
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
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
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
sdf_wavy_2d.w
int w
Definition: sdf_wavy_2d.py:6
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
MetaSpringBC::OpSpringALEMaterialLhs_dX_dx
LHS-operator for the pressure element (material configuration)
Definition: SpringElement.hpp:515
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
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
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
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::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
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
MetaSpringBC::OpSpringKs
LHS-operator for the springs element.
Definition: SpringElement.hpp:200
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::OpSpringALEMaterialLhs_dX_dX
LHS-operator for the pressure element (material configuration)
Definition: SpringElement.hpp:544
Generated by Doxygen 1.8.17 and hosted at University of Glasgow