v0.14.0
CellForces.hpp
Go to the documentation of this file.
1 /* This file is part of MoFEM.
2 
3  * MoFEM is free software: you can redistribute it and/or modify it under
4  * the terms of the GNU Lesser General Public License as published by the
5  * Free Software Foundation, either version 3 of the License, or (at your
6  * option) any later version.
7  *
8  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
9  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
10  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
11  * License for more details.
12  *
13  * You should have received a copy of the GNU Lesser General Public
14  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
15 
16 #ifndef __CELL_FORCES_HPP__
17 #define __CELL_FORCES_HPP__
18 
19 namespace CellEngineering {
20 
21 struct FatPrism : public MoFEM::FatPrismElementForcesAndSourcesCore {
22  FatPrism(MoFEM::Interface &m_field)
23  : MoFEM::FatPrismElementForcesAndSourcesCore(m_field) {}
24  int getRuleTrianglesOnly(int order) { return 2 * (order - 1); }
25  int getRuleThroughThickness(int order) { return 2 * order; }
26 };
27 
28 /**
29  * Common data structure used to pass data between operators
30  */
31 struct CommonData {
32 
33  MatrixDouble forceAtGaussPts;
34  VectorDouble forcePotentialAtGaussPoints;
35  MatrixDouble gradForceAtGaussPtrs;
36 
37  boost::shared_ptr<VectorDouble> dispX;
38  boost::shared_ptr<VectorDouble> dispY;
39 
40  CommonData()
41  : dispX(boost::shared_ptr<VectorDouble>(new VectorDouble())),
42  dispY(boost::shared_ptr<VectorDouble>(new VectorDouble())) {}
43 };
44 
45 /**
46  * Face element used to run operators
47  */
48 struct FaceElement : public MoFEM::FaceElementForcesAndSourcesCore {
49  FaceElement(MoFEM::Interface &m_field)
50  : MoFEM::FaceElementForcesAndSourcesCore(m_field) {}
51  int getRule(int order) { return 2 * order; };
52 };
53 
54 /**
55  * Operator to get measured displacements in X-direction on the face
56  */
57 struct OpGetDispX : public MoFEM::OpCalculateScalarFieldValues {
58  OpGetDispX(CommonData &common_data)
59  : MoFEM::OpCalculateScalarFieldValues("DISP_X", common_data.dispX) {}
60 };
61 
62 /**
63  * Operator to get measured displacements in Y-direction on the face
64  */
65 struct OpGetDispY : public MoFEM::OpCalculateScalarFieldValues {
66  OpGetDispY(CommonData &common_data)
67  : MoFEM::OpCalculateScalarFieldValues("DISP_Y", common_data.dispY) {}
68 };
69 
70 /**
71  * \brief Calculate and assemble D matrix
72  */
73 struct OpCellPotentialD
74  : public MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator {
75 
76  Mat A;
77  const double epsRho;
78 
79  OpCellPotentialD(Mat a, const double rho)
80  : MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator(
81  "RHO", "RHO", UserDataOperator::OPROWCOL),
82  A(a), epsRho(rho) {
83  sYmm = false;
84  }
85 
86  MatrixDouble C;
87 
88  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
89  EntityType col_type,
90  DataForcesAndSourcesCore::EntData &row_data,
91  DataForcesAndSourcesCore::EntData &col_data) {
92 
93  MoFEMFunctionBegin;
94 
95  int nb_row = row_data.getIndices().size();
96  int nb_col = col_data.getIndices().size();
97 
98  if (nb_row == 0)
99  MoFEMFunctionReturnHot(0);
100  if (nb_col == 0)
101  MoFEMFunctionReturnHot(0);
102 
103  const VectorInt &row_indices = row_data.getIndices();
104  const VectorInt &col_indices = col_data.getIndices();
105 
106  const int nb_gauss_pts = row_data.getN().size1();
107 
108  C.resize(nb_row, nb_col, false);
109  C.clear();
110  for (int gg = 0; gg != nb_gauss_pts; gg++) {
111  // get integration weights
112  const double val = getGaussPts()(2, gg) * getArea() * epsRho;
113  for (int rr = 0; rr != nb_row; rr++) {
114  const double diff_nr_x = row_data.getDiffN()(gg, 2 * rr + 0);
115  const double diff_nr_y = row_data.getDiffN()(gg, 2 * rr + 1);
116  for (int cc = 0; cc != nb_col; cc++) {
117  const double diff_nc_x = col_data.getDiffN()(gg, 2 * cc + 0);
118  const double diff_nc_y = col_data.getDiffN()(gg, 2 * cc + 1);
119  C(rr, cc) += val * (diff_nr_x * diff_nc_x + diff_nr_y * diff_nc_y);
120  }
121  }
122  }
123 
124  CHKERR MatSetValues(A, nb_row, &*row_indices.begin(), nb_col,
125  &*col_indices.begin(), &*C.data().begin(), ADD_VALUES);
126 
127  MoFEMFunctionReturn(0);
128  }
129 };
130 
131 /**
132  * \brief Calculate and assemble B matrix
133 
134  \f[
135  \mathbf{B} = \int_{S_0}
136  \mathbf{N}^\textrm{T}_\textrm{u} \frac{\partial \mathbf{N}_\rho}{\partial
137  \mathbf{x}} \textrm{d}S \f]
138 
139  */
140 struct OpCellPotentialB
141  : public MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator {
142 
143  Mat A;
144 
145  OpCellPotentialB(Mat a, const string field_name)
146  : MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator(
147  "RHO", field_name, UserDataOperator::OPROWCOL),
148  A(a) {
149  sYmm = false;
150  }
151 
152  MatrixDouble C;
153 
154  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
155  EntityType col_type,
156  DataForcesAndSourcesCore::EntData &row_data,
157  DataForcesAndSourcesCore::EntData &col_data) {
158 
159  MoFEMFunctionBegin;
160 
161  int nb_row = row_data.getIndices().size();
162  int nb_col = col_data.getIndices().size();
163 
164  if (nb_row == 0)
165  MoFEMFunctionReturnHot(0);
166  if (nb_col == 0)
167  MoFEMFunctionReturnHot(0);
168 
169  const VectorInt &row_indices = row_data.getIndices();
170  const VectorInt &col_indices = col_data.getIndices();
171 
172  const int nb_gauss_pts = row_data.getN().size1();
173 
174  C.resize(nb_row, nb_col, false);
175  C.clear();
176  for (int gg = 0; gg != nb_gauss_pts; gg++) {
177  // get integration weights
178  const double val = getGaussPts()(2, gg) * getArea();
179  for (int rr = 0; rr != nb_row; rr++) {
180  const double diff_nr_x = row_data.getDiffN()(gg, 2 * rr + 0);
181  const double diff_nr_y = row_data.getDiffN()(gg, 2 * rr + 1);
182  for (int cc = 0; cc != nb_col / 3; cc++) {
183  const double nc = col_data.getN(gg)[cc];
184  C(rr, 3 * cc + 0) -= val * diff_nr_x * nc;
185  C(rr, 3 * cc + 1) -= val * diff_nr_y * nc;
186  }
187  }
188  }
189 
190  CHKERR MatSetValues(A, nb_row, &*row_indices.begin(), nb_col,
191  &*col_indices.begin(), &*C.data().begin(), ADD_VALUES);
192 
193  MoFEMFunctionReturn(0);
194  }
195 };
196 
197 /**
198  * \brief Calculate and assemble S matrix
199 
200  */
201 struct OpCellS
202  : public MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator {
203 
204  Mat A;
205  const double epsU;
206 
207  OpCellS(Mat a, const double eps_u)
208  : MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator(
209  "UPSILON", "U", UserDataOperator::OPROWCOL),
210  A(a), epsU(eps_u) {
211  sYmm = false;
212  }
213 
214  MatrixDouble C;
215 
216  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
217  EntityType col_type,
218  DataForcesAndSourcesCore::EntData &row_data,
219  DataForcesAndSourcesCore::EntData &col_data) {
220 
221  MoFEMFunctionBegin;
222 
223  int nb_row = row_data.getIndices().size();
224  int nb_col = col_data.getIndices().size();
225 
226  if (nb_row == 0)
227  MoFEMFunctionReturnHot(0);
228  if (nb_col == 0)
229  MoFEMFunctionReturnHot(0);
230 
231  const VectorInt &row_indices = row_data.getIndices();
232  const VectorInt &col_indices = col_data.getIndices();
233 
234  const int nb_gauss_pts = row_data.getN().size1();
235 
236  C.resize(nb_row, nb_col, false);
237  C.clear();
238  for (int gg = 0; gg != nb_gauss_pts; gg++) {
239  // get integration weights
240  const double val = getGaussPts()(2, gg) * getArea() * pow(epsU, -1);
241  for (int rr = 0; rr != nb_row / 3; rr++) {
242  const double nr = row_data.getN(gg)[rr];
243  for (int cc = 0; cc != nb_col / 3; cc++) {
244  const double nc = col_data.getN(gg)[cc];
245  const double c = val * nr * nc;
246  C(3 * rr + 0, 3 * cc + 0) += c;
247  C(3 * rr + 1, 3 * cc + 1) += c;
248  }
249  }
250  }
251 
252  CHKERR MatSetValues(A, nb_row, &*row_indices.begin(), nb_col,
253  &*col_indices.begin(), &*C.data().begin(), ADD_VALUES);
254 
255  MoFEMFunctionReturn(0);
256  }
257 };
258 
259 /**
260  * \brief Calculate and assemble g vector
261 
262  \f[
263  \mathbf{g} = \int_{S_1} \mathbf{N}^\textrm{T}_\textrm{u} \overline{\mathbf{u}}
264  \textrm{d}S, \f]
265 
266  */
267 struct OpCell_g
268  : public MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator {
269 
270  Vec F;
271  const double epsU;
272  CommonData &commonData;
273 
274  OpCell_g(Vec f, const double eps_u, CommonData &common_data)
275  : MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator(
276  "UPSILON", UserDataOperator::OPROW),
277  F(f), epsU(eps_u), commonData(common_data) {}
278 
279  VectorDouble nF;
280 
281  MoFEMErrorCode doWork(int side, EntityType type,
282  DataForcesAndSourcesCore::EntData &data) {
283 
284  MoFEMFunctionBegin;
285 
286  int nb_row = data.getIndices().size();
287  if (nb_row == 0)
288  MoFEMFunctionReturnHot(0);
289 
290  const VectorDouble &disp_x = *commonData.dispX;
291  const VectorDouble &disp_y = *commonData.dispY;
292 
293  const int nb_gauss_pts = data.getN().size1();
294 
295  nF.resize(nb_row, false);
296  nF.clear();
297  for (int gg = 0; gg != nb_gauss_pts; gg++) {
298  const double val = getGaussPts()(2, gg) * getArea() * pow(epsU, -1);
299  for (int rr = 0; rr != nb_row / 3; rr++) {
300  const double nr = val * data.getN(gg)[rr];
301  nF[3 * rr + 0] += nr * disp_x[gg];
302  nF[3 * rr + 1] += nr * disp_y[gg];
303  }
304  }
305  // assemble local face vector
306  CHKERR VecSetValues(F, nb_row, &data.getIndices()[0], &nF[0], ADD_VALUES);
307 
308  MoFEMFunctionReturn(0);
309  }
310 };
311 
312 /**
313  * \brief Calculate and assemble Z matrix
314  */
315 struct OpCellCurlD
316  : public MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator {
317 
318  Mat A;
319  const double epsRho;
320  const double epsL;
321 
322  OpCellCurlD(Mat a, const double eps_rho, const double eps_l)
323  : MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator(
324  "RHO", "RHO", UserDataOperator::OPROWCOL),
325  A(a), epsRho(eps_rho), epsL(eps_l) {
326  sYmm = false;
327  }
328 
329  MatrixDouble C;
330 
331  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
332  EntityType col_type,
333  DataForcesAndSourcesCore::EntData &row_data,
334  DataForcesAndSourcesCore::EntData &col_data) {
335 
336  MoFEMFunctionBegin;
337 
338  int nb_row = row_data.getIndices().size();
339  int nb_col = col_data.getIndices().size();
340 
341  if (nb_row == 0)
342  MoFEMFunctionReturnHot(0);
343  if (nb_col == 0)
344  MoFEMFunctionReturnHot(0);
345 
346  const VectorInt &row_indices = row_data.getIndices();
347  const VectorInt &col_indices = col_data.getIndices();
348 
349  const int nb_gauss_pts = row_data.getN().size1();
350 
351  auto t_diff_base_row = row_data.getFTensor2DiffN<3, 2>();
352  auto t_base_row = row_data.getFTensor1N<3>();
353 
354  FTensor::Index<'i', 3> i;
355  FTensor::Index<'j', 2> j;
356  FTensor::Number<0> N0;
357  FTensor::Number<1> N1;
358 
359  C.resize(nb_row, nb_col, false);
360  C.clear();
361  for (int gg = 0; gg != nb_gauss_pts; gg++) {
362  // get integration weights
363  const double val = getGaussPts()(2, gg) * getArea();
364  for (int rr = 0; rr != nb_row; rr++) {
365  const double curl_row = t_diff_base_row(1, 0) - t_diff_base_row(0, 1);
366  auto t_diff_base_col = col_data.getFTensor2DiffN<3, 2>(gg, 0);
367  auto t_base_col = col_data.getFTensor1N<3>(gg, 0);
368  for (int cc = 0; cc != nb_col; cc++) {
369  const double curl_col = t_diff_base_col(1, 0) - t_diff_base_col(0, 1);
370  const double diag = val * epsRho * t_base_row(i) * t_base_col(i);
371  const double curl =
372  epsL > 0 ? val * pow(epsL, -1) * curl_row * curl_col : 0;
373  C(rr, cc) += diag + curl;
374  ++t_diff_base_col;
375  ++t_base_col;
376  }
377  ++t_diff_base_row;
378  ++t_base_row;
379  }
380  }
381 
382  CHKERR MatSetValues(A, nb_row, &*row_indices.begin(), nb_col,
383  &*col_indices.begin(), &*C.data().begin(), ADD_VALUES);
384 
385  MoFEMFunctionReturn(0);
386  }
387 };
388 
389 /**
390  * \brief Calculate and assemble Z matrix
391  */
392 struct OpCellCurlB
393  : public MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator {
394 
395  Mat A;
396 
397  OpCellCurlB(Mat a, const string field_name)
398  : MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator(
399  "RHO", field_name, UserDataOperator::OPROWCOL),
400  A(a) {
401  sYmm = false;
402  }
403 
404  MatrixDouble C;
405 
406  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
407  EntityType col_type,
408  DataForcesAndSourcesCore::EntData &row_data,
409  DataForcesAndSourcesCore::EntData &col_data) {
410 
411  MoFEMFunctionBegin;
412 
413  int nb_row = row_data.getIndices().size();
414  int nb_col = col_data.getIndices().size();
415 
416  if (nb_row == 0)
417  MoFEMFunctionReturnHot(0);
418  if (nb_col == 0)
419  MoFEMFunctionReturnHot(0);
420 
421  const VectorInt &row_indices = row_data.getIndices();
422  const VectorInt &col_indices = col_data.getIndices();
423 
424  const int nb_gauss_pts = row_data.getN().size1();
425 
426  auto t_base_row = row_data.getFTensor1N<3>();
427 
428  C.resize(nb_row, nb_col, false);
429  C.clear();
430  for (int gg = 0; gg != nb_gauss_pts; gg++) {
431  // get integration weights
432  const double val = getGaussPts()(2, gg) * getArea();
433  for (int rr = 0; rr != nb_row; rr++) {
434  // This is scalar field, diff yield tensor1
435  auto t_base_col = col_data.getFTensor0N(gg, 0);
436  for (int cc = 0; cc != nb_col / 3; cc++) {
437  C(rr, 3 * cc + 0) -= val * t_base_row(0) * t_base_col;
438  C(rr, 3 * cc + 1) -= val * t_base_row(1) * t_base_col;
439  ++t_base_col;
440  }
441  ++t_base_row;
442  }
443  }
444 
445  CHKERR MatSetValues(A, nb_row, &*row_indices.begin(), nb_col,
446  &*col_indices.begin(), &*C.data().begin(), ADD_VALUES);
447 
448  MoFEMFunctionReturn(0);
449  }
450 };
451 
452 /**
453  * \brief Post-process tractions
454  *
455  * \todo Should be done with Prisms, no need to have three kinds of elements,
456  * tets for solid, prisms for constrains and triangles to do post-processing. We
457  * can do post-pocessing with prisms.
458  *
459  */
460 struct OpVirtualPotentialRho
461  : public MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator {
462 
463  CommonData &commonData;
464 
465  OpVirtualPotentialRho(std::string field_name, CommonData &common_data)
466  : MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator(
467  field_name, UserDataOperator::OPROW),
468  commonData(common_data) {}
469 
470  /**
471  * Real work is done here, i.e. assemble tractions
472  */
473  MoFEMErrorCode doWork(int side, EntityType type,
474  DataForcesAndSourcesCore::EntData &data) {
475  //
476  //
477  MoFEMFunctionBegin;
478 
479  if (data.getFieldData().empty())
480  MoFEMFunctionReturnHot(0);
481 
482  const int nb_gauss_pts = data.getN().size1();
483  const int nb_dofs = data.getFieldData().size();
484 
485  VectorDouble *vec;
486  MatrixDouble *mat;
487  MatrixDouble *grad;
488 
489  if (rowFieldName == "RHO") {
490  vec = &(commonData.forcePotentialAtGaussPoints);
491  mat = &commonData.forceAtGaussPts;
492  grad = &commonData.gradForceAtGaussPtrs;
493  }
494 
495  if (type == MBVERTEX) {
496  mat->resize(nb_gauss_pts, 2, false);
497  mat->clear();
498  vec->resize(nb_gauss_pts, false);
499  vec->clear();
500  grad->resize(0, 0, false);
501  }
502 
503  // calculate forces at integration points
504  for (unsigned int gg = 0; gg < nb_gauss_pts; gg++) {
505  for (int rr = 0; rr != nb_dofs; rr++) {
506  const double diff_base_x = data.getDiffN()(gg, 2 * rr + 0);
507  const double diff_base_y = data.getDiffN()(gg, 2 * rr + 1);
508  const double val = data.getFieldData()[rr];
509  (*mat)(gg, 0) += diff_base_x * val;
510  (*mat)(gg, 1) += diff_base_y * val;
511  const double base = data.getN()(gg, rr);
512  (*vec)[gg] += base * val;
513  }
514  }
515 
516  MoFEMFunctionReturn(0);
517  }
518 };
519 
520 /**
521  * \brief Post-process tractions
522  *
523  * \todo Should be done with Prisms, no need to have three kinds of elements,
524  * tets for solid, prisms for constrains and triangles to do post-processing. We
525  * can do post-pocessing with prisms.
526  *
527  */
528 struct OpVirtualCurlRho
529  : public MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator {
530 
531  CommonData &commonData;
532 
533  OpVirtualCurlRho(std::string field_name, CommonData &common_data)
534  : MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator(
535  field_name, UserDataOperator::OPROW),
536  commonData(common_data) {}
537 
538  /**
539  * Real work is done here, i.e. assemble tractions
540  */
541  MoFEMErrorCode doWork(int side, EntityType type,
542  DataForcesAndSourcesCore::EntData &data) {
543  //
544  //
545  MoFEMFunctionBegin;
546 
547  if (data.getFieldData().empty())
548  MoFEMFunctionReturnHot(0);
549 
550  const int nb_gauss_pts = data.getN().size1();
551  const int nb_dofs = data.getFieldData().size();
552 
553  VectorDouble *vec;
554  MatrixDouble *mat;
555  MatrixDouble *grad;
556 
557  if (rowFieldName == "RHO") {
558  mat = &commonData.forceAtGaussPts;
559  grad = &commonData.gradForceAtGaussPtrs;
560  vec = &commonData.forcePotentialAtGaussPoints;
561  }
562 
563  if (type == MBEDGE && side == 0) {
564  mat->resize(nb_gauss_pts, 2, false);
565  mat->clear();
566  grad->resize(nb_gauss_pts, 4, false);
567  grad->clear();
568  vec->resize(0, false);
569  }
570 
571  // calculate forces at integration points
572  auto t_base = data.getFTensor1N<3>();
573  auto t_diff_base_fun = data.getFTensor2DiffN<3, 2>();
574 
575  for (unsigned int gg = 0; gg < nb_gauss_pts; gg++) {
576  auto t_data = data.getFTensor0FieldData();
577  for (int rr = 0; rr != nb_dofs; rr++) {
578  const double base_x = t_base(0);
579  const double base_y = t_base(1);
580  (*mat)(gg, 0) += base_x * t_data;
581  (*mat)(gg, 1) += base_y * t_data;
582  for (int i = 0; i != 2; i++) {
583  for (int j = 0; j != 2; j++) {
584  (*grad)(gg, 2 * i + j) += t_diff_base_fun(i, j) * t_data;
585  }
586  }
587  ++t_data;
588  ++t_base;
589  ++t_diff_base_fun;
590  }
591  }
592 
593  MoFEMFunctionReturn(0);
594  }
595 };
596 
597 /**
598  * \brief Shave results on mesh tags for post-processing
599  */
600 struct PostProcTraction
601  : public MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator {
602 
603  MoFEM::Interface &mField;
604  moab::Interface &postProcMesh;
605  std::vector<EntityHandle> &mapGaussPts;
606  CommonData &commonData;
607 
608  /**
609  * Constructor
610  */
611  PostProcTraction(MoFEM::Interface &m_field, moab::Interface &post_proc_mesh,
612  std::vector<EntityHandle> &map_gauss_pts,
613  CommonData &common_data)
614  : MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator(
615  "RHO", ForcesAndSourcesCore::UserDataOperator::OPROW),
616  mField(m_field), postProcMesh(post_proc_mesh),
617  mapGaussPts(map_gauss_pts), commonData(common_data) {}
618 
619  /**
620  * \brief Here real work is done
621  */
622  MoFEMErrorCode doWork(int side, EntityType type,
623  DataForcesAndSourcesCore::EntData &data) {
624  MoFEMFunctionBegin;
625 
626  if (type != MBVERTEX)
627  MoFEMFunctionReturnHot(0);
628 
629  double def_val[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
630  Tag th_force, th_force_potential, th_grad_force;
631  CHKERR postProcMesh.tag_get_handle("FORCE", 3, MB_TYPE_DOUBLE, th_force,
632  MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
633 
634  if (commonData.forcePotentialAtGaussPoints.size() > 0) {
635  CHKERR postProcMesh.tag_get_handle("FORCE_POTENTIAL", 1, MB_TYPE_DOUBLE,
636  th_force_potential,
637  MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
638  }
639  if (commonData.gradForceAtGaussPtrs.size1() > 0) {
640  CHKERR postProcMesh.tag_get_handle("FORCE_GRADIENT", 9, MB_TYPE_DOUBLE,
641  th_grad_force,
642  MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
643  }
644 
645  int nb_gauss_pts = data.getN().size1();
646  if (mapGaussPts.size() != (unsigned int)nb_gauss_pts) {
647  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
648  }
649 
650  for (int gg = 0; gg < nb_gauss_pts; gg++) {
651  const double force[] = {commonData.forceAtGaussPts(gg, 0),
652  commonData.forceAtGaussPts(gg, 1), 0};
653  CHKERR postProcMesh.tag_set_data(th_force, &mapGaussPts[gg], 1, force);
654  if (commonData.forcePotentialAtGaussPoints.size() > 0) {
655  CHKERR postProcMesh.tag_set_data(
656  th_force_potential, &mapGaussPts[gg], 1,
657  &commonData.forcePotentialAtGaussPoints[gg]);
658  }
659  if (commonData.gradForceAtGaussPtrs.size1() > 0) {
660  const double grad_force[] = {commonData.gradForceAtGaussPtrs(gg, 0),
661  commonData.gradForceAtGaussPtrs(gg, 1),
662  0,
663  commonData.gradForceAtGaussPtrs(gg, 2),
664  commonData.gradForceAtGaussPtrs(gg, 3),
665  0,
666  0,
667  0,
668  0};
669  CHKERR postProcMesh.tag_set_data(th_grad_force, &mapGaussPts[gg], 1,
670  grad_force);
671 
672  }
673  }
674 
675  MoFEMFunctionReturn(0);
676  }
677 };
678 
679 } // namespace CellEngineering
680 
681 #endif // __CELL_FORCES_HPP__
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
CellEngineering::OpCellCurlD::A
Mat A
Definition: CellForces.hpp:318
CellEngineering::OpCellCurlB::C
MatrixDouble C
Definition: CellForces.hpp:404
CellEngineering::PostProcTraction
Shave results on mesh tags for post-processing.
Definition: CellForces.hpp:600
CellEngineering
Definition: cell_forces.cpp:412
MoFEM::MatSetValues
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
Definition: EntitiesFieldData.hpp:1631
CellEngineering::CommonData::dispX
boost::shared_ptr< VectorDouble > dispX
Definition: CellForces.hpp:37
CellEngineering::OpCellS::C
MatrixDouble C
Definition: CellForces.hpp:214
CellEngineering::PostProcTraction::postProcMesh
moab::Interface & postProcMesh
Definition: CellForces.hpp:604
CellEngineering::OpVirtualPotentialRho::OpVirtualPotentialRho
OpVirtualPotentialRho(std::string field_name, CommonData &common_data)
Definition: CellForces.hpp:465
rho
double rho
Definition: plastic.cpp:191
CellEngineering::FaceElement
Definition: CellForces.hpp:48
CellEngineering::OpCellPotentialB::A
Mat A
Definition: CellForces.hpp:143
CellEngineering::OpVirtualCurlRho::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: CellForces.hpp:541
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
CellEngineering::CommonData::CommonData
CommonData()
Definition: CellForces.hpp:40
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
CellEngineering::OpCellCurlB
Calculate and assemble Z matrix.
Definition: CellForces.hpp:392
CellEngineering::OpCellS::epsU
const double epsU
Definition: CellForces.hpp:205
CellEngineering::FaceElement::getRule
int getRule(int order)
Definition: CellForces.hpp:51
MoFEM::ForcesAndSourcesCore::UserDataOperator::rowFieldName
std::string rowFieldName
Definition: ForcesAndSourcesCore.hpp:577
CellEngineering::OpVirtualCurlRho::OpVirtualCurlRho
OpVirtualCurlRho(std::string field_name, CommonData &common_data)
Definition: CellForces.hpp:533
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPROWCOL
@ OPROWCOL
operator doWork is executed on FE rows &columns
Definition: ForcesAndSourcesCore.hpp:569
CellEngineering::OpGetDispY::OpGetDispY
OpGetDispY(CommonData &common_data)
Definition: CellForces.hpp:66
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1576
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
CellEngineering::OpCellPotentialB::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Definition: CellForces.hpp:154
MoFEM::ForcesAndSourcesCore::UserDataOperator::getGaussPts
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Definition: ForcesAndSourcesCore.hpp:1235
MoFEM::ForcesAndSourcesCore::UserDataOperator::ForcesAndSourcesCore
friend class ForcesAndSourcesCore
Definition: ForcesAndSourcesCore.hpp:989
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
CellEngineering::OpVirtualCurlRho
Post-process tractions.
Definition: CellForces.hpp:528
CellEngineering::OpCellPotentialD::A
Mat A
Definition: CellForces.hpp:76
MoFEM::FatPrismElementForcesAndSourcesCore::FatPrismElementForcesAndSourcesCore
FatPrismElementForcesAndSourcesCore(Interface &m_field)
Definition: FatPrismElementForcesAndSourcesCore.cpp:11
CellEngineering::FatPrism::FatPrism
FatPrism(MoFEM::Interface &m_field)
Definition: CellForces.hpp:22
CellEngineering::OpCell_g::OpCell_g
OpCell_g(Vec f, const double eps_u, CommonData &common_data)
Definition: CellForces.hpp:274
FTensor::Number< 0 >
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
CellEngineering::OpCellCurlD::C
MatrixDouble C
Definition: CellForces.hpp:329
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
CellEngineering::OpCellCurlD::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Definition: CellForces.hpp:331
CellEngineering::OpCellPotentialD::epsRho
const double epsRho
Definition: CellForces.hpp:77
a
constexpr double a
Definition: approx_sphere.cpp:30
MoFEM::ForcesAndSourcesCore::UserDataOperator
Definition: ForcesAndSourcesCore.hpp:549
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
default operator for TRI element
Definition: FaceElementForcesAndSourcesCore.hpp:94
convert.type
type
Definition: convert.py:64
CellEngineering::OpCellCurlD::epsRho
const double epsRho
Definition: CellForces.hpp:319
CellEngineering::OpCellCurlD
Calculate and assemble Z matrix.
Definition: CellForces.hpp:315
CellEngineering::OpCellPotentialD
Calculate and assemble D matrix.
Definition: CellForces.hpp:73
CellEngineering::OpCellS::OpCellS
OpCellS(Mat a, const double eps_u)
Definition: CellForces.hpp:207
CellEngineering::OpCellS::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Definition: CellForces.hpp:216
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
CellEngineering::FaceElement::FaceElement
FaceElement(MoFEM::Interface &m_field)
Definition: CellForces.hpp:49
CellEngineering::OpVirtualPotentialRho::commonData
CommonData & commonData
Definition: CellForces.hpp:463
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
CellEngineering::OpVirtualPotentialRho::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: CellForces.hpp:473
CellEngineering::FatPrism
Definition: CellForces.hpp:21
CellEngineering::OpCell_g::commonData
CommonData & commonData
Definition: CellForces.hpp:272
CellEngineering::OpCellS
Calculate and assemble S matrix.
Definition: CellForces.hpp:201
MoFEM::FatPrismElementForcesAndSourcesCore
FatPrism finite element.
Definition: FatPrismElementForcesAndSourcesCore.hpp:31
CellEngineering::OpGetDispX
Definition: CellForces.hpp:57
CellEngineering::OpCellPotentialD::C
MatrixDouble C
Definition: CellForces.hpp:86
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
CellEngineering::OpVirtualCurlRho::commonData
CommonData & commonData
Definition: CellForces.hpp:531
CellEngineering::OpCell_g::nF
VectorDouble nF
Definition: CellForces.hpp:279
FaceElementForcesAndSourcesCore
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
FTensor::Index< 'i', 3 >
CellEngineering::PostProcTraction::mapGaussPts
std::vector< EntityHandle > & mapGaussPts
Definition: CellForces.hpp:605
CellEngineering::OpCellPotentialB::C
MatrixDouble C
Definition: CellForces.hpp:152
CellEngineering::CommonData::gradForceAtGaussPtrs
MatrixDouble gradForceAtGaussPtrs
Definition: CellForces.hpp:35
CellEngineering::FatPrism::getRuleThroughThickness
int getRuleThroughThickness(int order)
Definition: CellForces.hpp:25
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
MoFEM::Types::VectorInt
UBlasVector< int > VectorInt
Definition: Types.hpp:67
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
CellEngineering::OpCellPotentialB
Calculate and assemble B matrix.
Definition: CellForces.hpp:140
CellEngineering::OpCellCurlB::OpCellCurlB
OpCellCurlB(Mat a, const string field_name)
Definition: CellForces.hpp:397
CellEngineering::OpCellPotentialD::OpCellPotentialD
OpCellPotentialD(Mat a, const double rho)
Definition: CellForces.hpp:79
CellEngineering::OpGetDispX::OpGetDispX
OpGetDispX(CommonData &common_data)
Definition: CellForces.hpp:58
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
CellEngineering::CommonData
Definition: CellForces.hpp:31
CellEngineering::CommonData::forcePotentialAtGaussPoints
VectorDouble forcePotentialAtGaussPoints
Definition: CellForces.hpp:34
CellEngineering::OpGetDispY
Definition: CellForces.hpp:65
CellEngineering::OpCellCurlB::A
Mat A
Definition: CellForces.hpp:395
CellEngineering::OpCellPotentialD::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Definition: CellForces.hpp:88
CellEngineering::CommonData::forceAtGaussPts
MatrixDouble forceAtGaussPts
Definition: CellForces.hpp:33
CellEngineering::FatPrism::getRuleTrianglesOnly
int getRuleTrianglesOnly(int order)
Definition: CellForces.hpp:24
CellEngineering::OpCellPotentialB::OpCellPotentialB
OpCellPotentialB(Mat a, const string field_name)
Definition: CellForces.hpp:145
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
CellEngineering::OpCellCurlD::OpCellCurlD
OpCellCurlD(Mat a, const double eps_rho, const double eps_l)
Definition: CellForces.hpp:322
CellEngineering::PostProcTraction::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Here real work is done.
Definition: CellForces.hpp:622
CellEngineering::OpCell_g::epsU
const double epsU
Definition: CellForces.hpp:271
CellEngineering::OpCell_g::F
Vec F
Definition: CellForces.hpp:270
CellEngineering::OpCell_g::doWork
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: CellForces.hpp:281
CellEngineering::PostProcTraction::PostProcTraction
PostProcTraction(MoFEM::Interface &m_field, moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, CommonData &common_data)
Definition: CellForces.hpp:611
CellEngineering::PostProcTraction::mField
MoFEM::Interface & mField
Definition: CellForces.hpp:603
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator::getArea
double getArea()
get area of face
Definition: FaceElementForcesAndSourcesCore.hpp:239
CellEngineering::OpCell_g
Calculate and assemble g vector.
Definition: CellForces.hpp:267
CellEngineering::PostProcTraction::commonData
CommonData & commonData
Definition: CellForces.hpp:606
CellEngineering::OpVirtualPotentialRho
Post-process tractions.
Definition: CellForces.hpp:460
CellEngineering::CommonData::dispY
boost::shared_ptr< VectorDouble > dispY
Definition: CellForces.hpp:38
MoFEM::DataOperator::sYmm
bool sYmm
If true assume that matrix is symmetric structure.
Definition: DataOperators.hpp:82
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
CellEngineering::OpCellS::A
Mat A
Definition: CellForces.hpp:204
CellEngineering::OpCellCurlD::epsL
const double epsL
Definition: CellForces.hpp:320
CellEngineering::OpCellCurlB::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Definition: CellForces.hpp:406
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPROW
@ OPROW
operator doWork function is executed on FE rows
Definition: ForcesAndSourcesCore.hpp:567
Generated by Doxygen 1.8.17 and hosted at University of Glasgow