v0.9.1
PostProcOnRefMesh.cpp
Go to the documentation of this file.
1 /** \file PostProcOnRefMesh.cpp
2  * \brief Postprocess fields on refined mesh made for 10 Node tets
3  *
4  * Create refined mesh, without enforcing continuity between element. Calculate
5  * field values on nodes of that mesh.
6  * \ingroup mofem_fs_post_proc
7  */
8 
9 /*
10  * This file is part of MoFEM.
11  * MoFEM is free software: you can redistribute it and/or modify it under
12  * the terms of the GNU Lesser General Public License as published by the
13  * Free Software Foundation, either version 3 of the License, or (at your
14  * option) any later version.
15  *
16  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
17  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
19  * License for more details.
20  *
21  * You should have received a copy of the GNU Lesser General Public
22  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
23 
24 #include <MoFEM.hpp>
25 using namespace MoFEM;
26 using namespace boost::numeric;
27 #include <PostProcOnRefMesh.hpp>
28 
30  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
32 
33  if (data.getFieldData().size() == 0)
35 
36  if (V) {
37  vAlues.resize(data.getFieldData().size());
38  double *a;
39  CHKERR VecGetArray(V, &a);
40  VectorDofs::iterator it, hi_it;
41  it = data.getFieldDofs().begin();
42  hi_it = data.getFieldDofs().end();
43  for (int ii = 0; it != hi_it; it++, ii++) {
44  int local_idx = getFEMethod()
45  ->rowPtr->find((*it)->getGlobalUniqueId())
46  ->get()
47  ->getPetscLocalDofIdx();
48  vAlues[ii] = a[local_idx];
49  }
50  CHKERR VecRestoreArray(V, &a);
51  vAluesPtr = &vAlues;
52  } else {
53  vAluesPtr = &data.getFieldData();
54  }
55 
56  const MoFEM::FEDofEntity *dof_ptr = data.getFieldDofs()[0].get();
57  int rank = dof_ptr->getNbOfCoeffs();
58 
59  int tag_length = rank;
60  FieldSpace space = dof_ptr->getSpace();
61  switch (space) {
62  case L2:
63  case H1:
64  break;
65  case HCURL:
66  case HDIV:
67  tag_length *= 3;
68  break;
69  default:
70  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
71  "field with that space is not implemented");
72  }
73 
74  if (tag_length > 1 && tag_length < 3)
75  tag_length = 3;
76  else if (tag_length > 3 && tag_length < 9)
77  tag_length = 9;
78 
79  double def_VAL[tag_length];
80  bzero(def_VAL, tag_length * sizeof(double));
81  Tag th;
82  CHKERR postProcMesh.tag_get_handle(tagName.c_str(), tag_length,
83  MB_TYPE_DOUBLE, th,
84  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
85 
86  // zero tags, this for Vertex if H1 and TRI if Hdiv, EDGE for Hcurl
87  // no need for L2
88  const void *tags_ptr[mapGaussPts.size()];
89  int nb_gauss_pts = data.getN().size1();
90  if (mapGaussPts.size() != (unsigned int)nb_gauss_pts) {
91  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
92  "data inconsistency %d!=%d", mapGaussPts.size(), nb_gauss_pts);
93  }
94 
95  switch (space) {
96  case H1:
97  commonData.fieldMap[rowFieldName].resize(nb_gauss_pts);
98  if (type == MBVERTEX) {
99  for (int gg = 0; gg < nb_gauss_pts; gg++) {
100  rval = postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1, def_VAL);
101  (commonData.fieldMap[rowFieldName])[gg].resize(rank);
102  (commonData.fieldMap[rowFieldName])[gg].clear();
103  }
104  }
105  CHKERR postProcMesh.tag_get_by_ptr(th, &mapGaussPts[0], mapGaussPts.size(),
106  tags_ptr);
107  for (int gg = 0; gg < nb_gauss_pts; gg++) {
108  for (int rr = 0; rr < rank; rr++) {
109  const double val =
110  cblas_ddot((vAluesPtr->size() / rank), &(data.getN(gg)[0]), 1,
111  &((*vAluesPtr)[rr]), rank);
112  (commonData.fieldMap[rowFieldName])[gg][rr] =
113  ((double *)tags_ptr[gg])[rr] += val;
114  }
115  }
116  break;
117  case L2:
118  commonData.fieldMap[rowFieldName].resize(nb_gauss_pts);
119  for (int gg = 0; gg < nb_gauss_pts; gg++) {
120  CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1, def_VAL);
121  (commonData.fieldMap[rowFieldName])[gg].resize(rank);
122  (commonData.fieldMap[rowFieldName])[gg].clear();
123  }
124  CHKERR postProcMesh.tag_get_by_ptr(th, &mapGaussPts[0], mapGaussPts.size(),
125  tags_ptr);
126  for (int gg = 0; gg < nb_gauss_pts; gg++) {
127  bzero((double *)tags_ptr[gg], sizeof(double) * tag_length);
128  for (int rr = 0; rr < rank; rr++) {
129  const double val =
130  cblas_ddot((vAluesPtr->size() / rank), &(data.getN(gg)[0]), 1,
131  &((*vAluesPtr)[rr]), rank);
132  (commonData.fieldMap[rowFieldName])[gg][rr] =
133  ((double *)tags_ptr[gg])[rr] = val;
134  }
135  }
136  break;
137  case HCURL:
138  // FIXME: fieldMap not set
139  if (type == MBEDGE && side == 0) {
140  for (int gg = 0; gg < nb_gauss_pts; gg++) {
141  CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1, def_VAL);
142  }
143  }
144  CHKERR postProcMesh.tag_get_by_ptr(th, &mapGaussPts[0], mapGaussPts.size(),
145  tags_ptr);
146  {
148  auto t_n_hcurl = data.getFTensor1N<3>();
149  for (int gg = 0; gg != nb_gauss_pts; gg++) {
150  double *ptr = &((double *)tags_ptr[gg])[0];
151  int ll = 0;
152  for (; ll != static_cast<int>(vAluesPtr->size() / rank); ll++) {
153  FTensor::Tensor1<double *, 3> t_tag_val(ptr, &ptr[1], &ptr[2], 3);
154  for (int rr = 0; rr != rank; rr++) {
155  const double dof_val = (*vAluesPtr)[ll * rank + rr];
156  t_tag_val(i) += dof_val * t_n_hcurl(i);
157  ++t_tag_val;
158  }
159  ++t_n_hcurl;
160  }
161  for (; ll != static_cast<int>(data.getN().size2() / 3); ll++) {
162  ++t_n_hcurl;
163  }
164  }
165  }
166  break;
167  case HDIV:
168  // FIXME: fieldMap not set
169  if (type == MBTRI && side == 0) {
170  for (int gg = 0; gg < nb_gauss_pts; gg++) {
171  CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1, def_VAL);
172  }
173  }
174  CHKERR postProcMesh.tag_get_by_ptr(th, &mapGaussPts[0], mapGaussPts.size(),
175  tags_ptr);
176  {
178  auto t_n_hdiv = data.getFTensor1N<3>();
179  for (int gg = 0; gg != nb_gauss_pts; gg++) {
180  double *ptr = &((double *)tags_ptr[gg])[0];
181  int ll = 0;
182  for (; ll != static_cast<int>(vAluesPtr->size() / rank); ll++) {
183  FTensor::Tensor1<double *, 3> t_tag_val(ptr, &ptr[1], &ptr[2], 3);
184  for (int rr = 0; rr != rank; rr++) {
185  const double dof_val = (*vAluesPtr)[ll * rank + rr];
186  t_tag_val(i) += dof_val * t_n_hdiv(i);
187  ++t_tag_val;
188  }
189  ++t_n_hdiv;
190  }
191  for (; ll != static_cast<int>(data.getN().size2() / 3); ll++) {
192  ++t_n_hdiv;
193  }
194  }
195  }
196  break;
197  default:
198  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
199  "field with that space is not implemented");
200  }
201 
203 }
204 
206  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
208 
209  if (data.getFieldData().size() == 0)
211  if (V) {
212  vAlues.resize(data.getFieldData().size());
213  double *a;
214  CHKERR VecGetArray(V, &a);
215  VectorDofs::iterator it, hi_it;
216  it = data.getFieldDofs().begin();
217  hi_it = data.getFieldDofs().end();
218  for (int ii = 0; it != hi_it; it++, ii++) {
219  int local_idx = getFEMethod()
220  ->rowPtr->find((*it)->getGlobalUniqueId())
221  ->get()
222  ->getPetscLocalDofIdx();
223  vAlues[ii] = a[local_idx];
224  }
225  CHKERR VecRestoreArray(V, &a);
226  vAluesPtr = &vAlues;
227  } else {
228  vAluesPtr = &data.getFieldData();
229  }
230 
231  const MoFEM::FEDofEntity *dof_ptr = data.getFieldDofs()[0].get();
232  int rank = dof_ptr->getNbOfCoeffs();
233 
234  int tag_length = rank * 3;
235  FieldSpace space = dof_ptr->getSpace();
236  switch (space) {
237  case L2:
238  case H1:
239  break;
240  case HCURL:
241  case HDIV:
242  tag_length *= 3;
243  break;
244  default:
245  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
246  "field with that space is not implemented");
247  }
248 
249  double def_VAL[tag_length];
250  bzero(def_VAL, tag_length * sizeof(double));
251  Tag th;
252  CHKERR postProcMesh.tag_get_handle(tagName.c_str(), tag_length,
253  MB_TYPE_DOUBLE, th,
254  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
255 
256  // zero tags, this for Vertex if H1 and TRI if Hdiv, EDGE for Hcurl
257  // no need for L2
258  const void *tags_ptr[mapGaussPts.size()];
259  int nb_gauss_pts = data.getN().size1();
260  if (mapGaussPts.size() != (unsigned int)nb_gauss_pts) {
261  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
262  }
263 
264  switch (space) {
265  case H1:
266  commonData.gradMap[rowFieldName].resize(nb_gauss_pts);
267  if (type == MBVERTEX) {
268  for (int gg = 0; gg < nb_gauss_pts; gg++) {
269  CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1, def_VAL);
270  (commonData.gradMap[rowFieldName])[gg].resize(rank, 3);
271  (commonData.gradMap[rowFieldName])[gg].clear();
272  }
273  }
274  CHKERR postProcMesh.tag_get_by_ptr(th, &mapGaussPts[0], mapGaussPts.size(),
275  tags_ptr);
276  for (int gg = 0; gg < nb_gauss_pts; gg++) {
277  for (int rr = 0; rr < rank; rr++) {
278  for (int dd = 0; dd < 3; dd++) {
279  for (unsigned int dof = 0; dof < (vAluesPtr->size() / rank); dof++) {
280  const double val =
281  data.getDiffN(gg)(dof, dd) * (*vAluesPtr)[rank * dof + rr];
282  (commonData.gradMap[rowFieldName])[gg](rr, dd) =
283  ((double *)tags_ptr[gg])[rank * rr + dd] += val;
284  }
285  }
286  }
287  }
288  break;
289  case L2:
290  commonData.gradMap[rowFieldName].resize(nb_gauss_pts);
291  for (int gg = 0; gg < nb_gauss_pts; gg++) {
292  CHKERR postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1, def_VAL);
293  (commonData.gradMap[rowFieldName])[gg].resize(rank, 3);
294  (commonData.gradMap[rowFieldName])[gg].clear();
295  }
296  CHKERR postProcMesh.tag_get_by_ptr(th, &mapGaussPts[0], mapGaussPts.size(),
297  tags_ptr);
298  for (int gg = 0; gg < nb_gauss_pts; gg++) {
299  for (int rr = 0; rr < rank; rr++) {
300  for (int dd = 0; dd < 3; dd++) {
301  for (unsigned int dof = 0; dof < (vAluesPtr->size() / rank); dof++) {
302  const double val =
303  data.getDiffN(gg)(dof, dd) * (*vAluesPtr)[rank * dof + rr];
304  (commonData.gradMap[rowFieldName])[gg](rr, dd) =
305  ((double *)tags_ptr[gg])[rank * rr + dd] += val;
306  }
307  }
308  }
309  }
310  break;
311  default:
312  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
313  "field with that space is not implemented");
314  }
315 
317 }
318 
322 }
323 
325  int order_triangles_only) {
327  // if(gaussPtsTrianglesOnly.size1()==0 || gaussPtsTrianglesOnly.size2()==0) {
328  // SETERRQ(PETSC_COMM_SELF,MOFEM_DATA_INCONSISTENCY,"post-process mesh not
329  // generated");
330  // }
331 
332  // FIXME: Refinement not implement and inefficient implementation, looks ugly.
333 
334  //
335  const EntityHandle *conn;
336  int num_nodes;
337  EntityHandle prism;
338  const int nb_on_triangle = 6;
339  const int nb_through_thickness = 3;
340 
341  if (elementsMap.find(numeredEntFiniteElementPtr->getEnt()) !=
342  elementsMap.end()) {
343  prism = elementsMap[numeredEntFiniteElementPtr->getEnt()];
344  } else {
345  MatrixDouble gauss_pts_triangles_only, gauss_pts_through_thickness;
346  {
347  gauss_pts_triangles_only.resize(3, 3, false);
348  gauss_pts_triangles_only.clear();
349  gauss_pts_triangles_only(0, 0) = 0;
350  gauss_pts_triangles_only(1, 0) = 0;
351  gauss_pts_triangles_only(0, 1) = 1;
352  gauss_pts_triangles_only(1, 1) = 0;
353  gauss_pts_triangles_only(0, 2) = 0;
354  gauss_pts_triangles_only(1, 2) = 1;
355  gauss_pts_through_thickness.resize(2, 2, false);
356  gauss_pts_through_thickness(0, 0) = 0;
357  gauss_pts_through_thickness(0, 1) = 1;
358  }
359  ublas::vector<EntityHandle> prism_conn(6);
360  VectorDouble coords(3);
361  for (int ggf = 0; ggf != 3; ggf++) {
362  double ksi = gauss_pts_triangles_only(0, ggf);
363  double eta = gauss_pts_triangles_only(1, ggf);
364  coords[0] = ksi;
365  coords[1] = eta;
366  for (int ggt = 0; ggt != 2; ggt++) {
367  double zeta = gauss_pts_through_thickness(0, ggt);
368  coords[2] = zeta;
369  int side = ggt * 3 + ggf;
370  CHKERR postProcMesh.create_vertex(&coords[0], prism_conn[side]);
371  }
372  }
373  CHKERR postProcMesh.create_element(MBPRISM, &prism_conn[0], 6, prism);
374 
375  elementsMap[numeredEntFiniteElementPtr->getEnt()] = prism;
376  // Range faces;
377  // CHKERR postProcMesh.get_adjacencies(&prism,1,2,true,faces);
378  Range edges;
379  CHKERR postProcMesh.get_adjacencies(&prism, 1, 1, true, edges);
380  EntityHandle meshset;
381  CHKERR postProcMesh.create_meshset(MESHSET_SET, meshset);
382  CHKERR postProcMesh.add_entities(meshset, &prism, 1);
383  // CHKERR postProcMesh.add_entities(meshset,faces);
384  CHKERR postProcMesh.add_entities(meshset, edges);
385  if (tenNodesPostProcTets) {
386  CHKERR postProcMesh.convert_entities(meshset, true, false, false);
387  }
388  CHKERR postProcMesh.delete_entities(&meshset, 1);
389  CHKERR postProcMesh.delete_entities(edges);
390  // CHKERR postProcMesh.delete_entities(faces);
391 
392  CHKERR mField.get_moab().get_connectivity(
393  numeredEntFiniteElementPtr->getEnt(), conn, num_nodes, true);
394  MatrixDouble coords_prism_global;
395  coords_prism_global.resize(num_nodes, 3, false);
396  CHKERR mField.get_moab().get_coords(conn, num_nodes,
397  &coords_prism_global(0, 0));
398 
399  CHKERR postProcMesh.get_connectivity(prism, conn, num_nodes, false);
400  MatrixDouble coords_prism_local;
401  coords_prism_local.resize(num_nodes, 3, false);
402  CHKERR postProcMesh.get_coords(conn, num_nodes, &coords_prism_local(0, 0));
403 
404  if (gaussPtsThroughThickness.size2() != nb_through_thickness) {
405  pointsMap.clear();
406  gaussPtsTrianglesOnly.resize(3, nb_on_triangle, false);
407  gaussPtsTrianglesOnly.clear();
408  gaussPtsThroughThickness.resize(2, nb_through_thickness, false);
409  gaussPtsThroughThickness.clear();
410  const double eps = 1e-6;
411  int ggf = 0, ggt = 0;
412  for (int nn = 0; nn != num_nodes; nn++) {
413  double ksi = coords_prism_local(nn, 0);
414  double eta = coords_prism_local(nn, 1);
415  double zeta = coords_prism_local(nn, 2);
416  pointsMap.insert(PointsMap3D(ksi * 100, eta * 100, zeta * 100, nn));
417  if (fabs(zeta) < eps) {
418  gaussPtsTrianglesOnly(0, ggf) = ksi;
419  gaussPtsTrianglesOnly(1, ggf) = eta;
420  ggf++;
421  }
422  if (fabs(ksi) < eps && fabs(eta) < eps) {
423  gaussPtsThroughThickness(0, ggt) = zeta;
424  ggt++;
425  }
426  }
427  }
428 
429  for (int nn = 0; nn != num_nodes; nn++) {
430  const double ksi = coords_prism_local(nn, 0);
431  const double eta = coords_prism_local(nn, 1);
432  const double zeta = coords_prism_local(nn, 2);
433  const double n0 = N_MBTRI0(ksi, eta);
434  const double n1 = N_MBTRI1(ksi, eta);
435  const double n2 = N_MBTRI2(ksi, eta);
436  const double e0 = N_MBEDGE0(zeta);
437  const double e1 = N_MBEDGE1(zeta);
438  double coords_global[3];
439  for (int dd = 0; dd != 3; dd++) {
440  coords_global[dd] = n0 * e0 * coords_prism_global(0, dd) +
441  n1 * e0 * coords_prism_global(1, dd) +
442  n2 * e0 * coords_prism_global(2, dd) +
443  n0 * e1 * coords_prism_global(3, dd) +
444  n1 * e1 * coords_prism_global(4, dd) +
445  n2 * e1 * coords_prism_global(5, dd);
446  }
447  // cerr << coords_global[0] << " " << coords_global[1] << " " <<
448  // coords_global[2] << endl;
449  CHKERR postProcMesh.set_coords(&conn[nn], 1, coords_global);
450  }
451  }
452 
453  mapGaussPts.clear();
454  mapGaussPts.resize(nb_through_thickness * nb_on_triangle);
455  CHKERR postProcMesh.get_connectivity(prism, conn, num_nodes, false);
456  fill(mapGaussPts.begin(), mapGaussPts.end(), 0);
457  {
458  int gg = 0;
459  for (unsigned int ggf = 0; ggf != gaussPtsTrianglesOnly.size2(); ggf++) {
460  const double ksi = gaussPtsTrianglesOnly(0, ggf);
461  const double eta = gaussPtsTrianglesOnly(1, ggf);
462  for (unsigned int ggt = 0; ggt != gaussPtsThroughThickness.size2();
463  ggt++, gg++) {
464  const double zeta = gaussPtsThroughThickness(0, ggt);
465  PointsMap3D_multiIndex::iterator it;
466  it =
467  pointsMap.find(boost::make_tuple(ksi * 100, eta * 100, zeta * 100));
468  if (it != pointsMap.end()) {
469  mapGaussPts[gg] = conn[it->nN];
470  }
471  }
472  }
473  }
474 
476 }
477 
479  int order_thickness) {
481  if (gaussPtsThroughThickness.size1() == 0 ||
482  gaussPtsThroughThickness.size2() == 0) {
483  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
484  "post-process mesh not generated");
485  }
487 }
488 
491  // MoAB
492  ParallelComm *pcomm_post_proc_mesh =
493  ParallelComm::get_pcomm(&postProcMesh, MYPCOMM_INDEX);
494  if (pcomm_post_proc_mesh != NULL) {
495  delete pcomm_post_proc_mesh;
496  }
497  // CHKERR postProcMesh.delete_mesh();
499 }
500 
503  ParallelComm *pcomm =
504  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
505  ParallelComm *pcomm_post_proc_mesh =
506  ParallelComm::get_pcomm(&postProcMesh, MYPCOMM_INDEX);
507  if (pcomm_post_proc_mesh == NULL) {
508  pcomm_post_proc_mesh = new ParallelComm(&postProcMesh, mField.get_comm());
509  }
510  Range prims;
511  CHKERR postProcMesh.get_entities_by_type(0, MBPRISM, prims, false);
512  // std::cerr << "total prims size " << prims.size() << std::endl;
513  int rank = pcomm->rank();
514  Range::iterator pit = prims.begin();
515  for (; pit != prims.end(); pit++) {
516  CHKERR postProcMesh.tag_set_data(pcomm_post_proc_mesh->part_tag(), &*pit, 1,
517  &rank);
518  }
519  CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
521 }
522 
525 
526  auto generate_tri = [&](auto &gauss_pts) {
528  gauss_pts.resize(3, 3, false);
529  gauss_pts.clear();
530  gauss_pts(0, 0) = 0;
531  gauss_pts(1, 0) = 0;
532  gauss_pts(0, 1) = 1;
533  gauss_pts(1, 1) = 0;
534  gauss_pts(0, 2) = 0;
535  gauss_pts(1, 2) = 1;
536 
537  moab::Core core_ref;
538  moab::Interface &moab_ref = core_ref;
539  const EntityHandle *conn;
540  int num_nodes;
541  EntityHandle tri_conn[3];
542  MatrixDouble coords(6, 3);
543  for (int gg = 0; gg != 3; gg++) {
544  coords(gg, 0) = gauss_pts(0, gg);
545  coords(gg, 1) = gauss_pts(1, gg);
546  coords(gg, 2) = 0;
547  CHKERR moab_ref.create_vertex(&coords(gg, 0), tri_conn[gg]);
548  }
549 
550  EntityHandle tri;
551  CHKERR moab_ref.create_element(MBTRI, tri_conn, 3, tri);
552 
553  if (sixNodePostProcTris) {
554  Range edges;
555  CHKERR moab_ref.get_adjacencies(&tri, 1, 1, true, edges);
556  EntityHandle meshset;
557  CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
558  CHKERR moab_ref.add_entities(meshset, &tri, 1);
559  CHKERR moab_ref.add_entities(meshset, edges);
560  CHKERR moab_ref.convert_entities(meshset, true, false, false);
561  CHKERR moab_ref.get_connectivity(tri, conn, num_nodes, false);
562  CHKERR moab_ref.get_coords(conn, num_nodes, &coords(0, 0));
563  gauss_pts.resize(3, num_nodes, false);
564  gauss_pts.clear();
565  for (int nn = 0; nn < 3; nn++) {
566  gauss_pts(0, nn) = coords(nn, 0);
567  gauss_pts(1, nn) = coords(nn, 1);
568  gauss_pts(0, 3 + nn) = coords(3 + nn, 0);
569  gauss_pts(1, 3 + nn) = coords(3 + nn, 1);
570  }
571  } else {
572  CHKERR moab_ref.get_connectivity(tri, conn, num_nodes, false);
573  CHKERR moab_ref.get_coords(conn, num_nodes, &coords(0, 0));
574  gauss_pts.resize(3, num_nodes, false);
575  gauss_pts.clear();
576  for (int nn = 0; nn < 3; nn++) {
577  gauss_pts(0, nn) = coords(nn, 0);
578  gauss_pts(1, nn) = coords(nn, 1);
579  }
580  }
582  };
583 
584  auto generate_quad = [&](auto &gauss_pts) {
586 
587  gauss_pts.resize(3, 4, false);
588  gauss_pts.clear();
589 
590  gauss_pts(0, 0) = 0;
591  gauss_pts(1, 0) = 0;
592  gauss_pts(0, 1) = 1;
593  gauss_pts(1, 1) = 0;
594  gauss_pts(0, 2) = 1;
595  gauss_pts(1, 2) = 1;
596  gauss_pts(0, 3) = 0;
597  gauss_pts(1, 3) = 1;
598 
599  moab::Core core_ref;
600  moab::Interface &moab_ref = core_ref;
601 
602  const EntityHandle *conn;
603  int num_nodes;
604  std::array<EntityHandle, 4> quad_conn;
605  MatrixDouble coords(8, 3);
606  for (int gg = 0; gg != 4; gg++) {
607  coords(gg, 0) = gauss_pts(0, gg);
608  coords(gg, 1) = gauss_pts(1, gg);
609  coords(gg, 2) = 0;
610  CHKERR moab_ref.create_vertex(&coords(gg, 0), quad_conn[gg]);
611  }
612 
613  EntityHandle quad;
614  CHKERR moab_ref.create_element(MBQUAD, quad_conn.data(), 4, quad);
615 
616  if (sixNodePostProcTris) {
617  Range edges;
618  CHKERR moab_ref.get_adjacencies(&quad, 1, 1, true, edges);
619  EntityHandle meshset;
620  CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
621  CHKERR moab_ref.add_entities(meshset, &quad, 1);
622  CHKERR moab_ref.add_entities(meshset, edges);
623  CHKERR moab_ref.convert_entities(meshset, true, false, false);
624  CHKERR moab_ref.get_connectivity(quad, conn, num_nodes, false);
625  CHKERR moab_ref.get_coords(conn, num_nodes, &coords(0, 0));
626  gauss_pts.resize(3, num_nodes, false);
627  gauss_pts.clear();
628  for (int nn = 0; nn != 4; nn++) {
629  gauss_pts(0, nn) = coords(nn, 0);
630  gauss_pts(1, nn) = coords(nn, 1);
631  gauss_pts(0, 4 + nn) = coords(4 + nn, 0);
632  gauss_pts(1, 4 + nn) = coords(4 + nn, 1);
633  }
634  } else {
635  CHKERR moab_ref.get_connectivity(quad, conn, num_nodes, false);
636  CHKERR moab_ref.get_coords(conn, num_nodes, &coords(0, 0));
637  gauss_pts.resize(3, num_nodes, false);
638  gauss_pts.clear();
639  for (int nn = 0; nn != 4; nn++) {
640  gauss_pts(0, nn) = coords(nn, 0);
641  gauss_pts(1, nn) = coords(nn, 1);
642  }
643  }
644 
646  };
647 
648  CHKERR generate_tri(gaussPtsTri);
649  CHKERR generate_quad(gaussPtsQuad);
650 
652 }
653 
656  if (gaussPtsTri.size1() == 0 && gaussPtsQuad.size1() == 0)
657  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
658  "post-process mesh not generated");
659 
660  auto create_tri = [&](auto &gauss_pts) {
661  EntityHandle tri;
662  std::array<EntityHandle, 3> tri_conn;
663  MatrixDouble3by3 coords_tri(3, 3);
664  VectorDouble3 coords(3);
665  CHKERR mField.get_moab().get_connectivity(
666  numeredEntFiniteElementPtr->getEnt(), conn, num_nodes, true);
667  CHKERR mField.get_moab().get_coords(conn, num_nodes, &coords_tri(0, 0));
668  for (int gg = 0; gg != 3; gg++) {
669  double ksi = gauss_pts(0, gg);
670  double eta = gauss_pts(1, gg);
671  double n0 = N_MBTRI0(ksi, eta);
672  double n1 = N_MBTRI1(ksi, eta);
673  double n2 = N_MBTRI2(ksi, eta);
674  double x =
675  n0 * coords_tri(0, 0) + n1 * coords_tri(1, 0) + n2 * coords_tri(2, 0);
676  double y =
677  n0 * coords_tri(0, 1) + n1 * coords_tri(1, 1) + n2 * coords_tri(2, 1);
678  double z =
679  n0 * coords_tri(0, 2) + n1 * coords_tri(1, 2) + n2 * coords_tri(2, 2);
680  coords[0] = x;
681  coords[1] = y;
682  coords[2] = z;
683  CHKERR postProcMesh.create_vertex(&coords[0], tri_conn[gg]);
684  }
685  CHKERR postProcMesh.create_element(MBTRI, &tri_conn[0], 3, tri);
686  return tri;
687  };
688 
689  auto create_quad = [&](auto &gauss_pts) {
690  EntityHandle quad;
691  std::array<EntityHandle, 4> quad_conn;
692  MatrixDouble coords_quad(4, 3);
693  VectorDouble coords(4);
694  CHKERR mField.get_moab().get_connectivity(
695  numeredEntFiniteElementPtr->getEnt(), conn, num_nodes, true);
696  CHKERR mField.get_moab().get_coords(conn, num_nodes, &coords_quad(0, 0));
697  for (int gg = 0; gg != 4; gg++) {
698  double ksi = gauss_pts(0, gg);
699  double eta = gauss_pts(1, gg);
700  double n0 = N_MBQUAD0(ksi, eta);
701  double n1 = N_MBQUAD1(ksi, eta);
702  double n2 = N_MBQUAD2(ksi, eta);
703  double n3 = N_MBQUAD3(ksi, eta);
704  double x = n0 * coords_quad(0, 0) + n1 * coords_quad(1, 0) +
705  n2 * coords_quad(2, 0) + n3 * coords_quad(3, 0);
706  double y = n0 * coords_quad(0, 1) + n1 * coords_quad(1, 1) +
707  n2 * coords_quad(2, 1) + n3 * coords_quad(3, 1);
708  double z = n0 * coords_quad(0, 2) + n1 * coords_quad(1, 2) +
709  n2 * coords_quad(2, 2) + n3 * coords_quad(3, 2);
710  coords[0] = x;
711  coords[1] = y;
712  coords[2] = z;
713  CHKERR postProcMesh.create_vertex(&coords[0], quad_conn[gg]);
714  }
715  CHKERR postProcMesh.create_element(MBQUAD, &quad_conn[0], 4, quad);
716  return quad;
717  };
718 
719  auto add_mid_nodes = [&](auto tri) {
721  Range edges;
722  CHKERR postProcMesh.get_adjacencies(&tri, 1, 1, true, edges);
723  EntityHandle meshset;
724  CHKERR postProcMesh.create_meshset(MESHSET_SET, meshset);
725  CHKERR postProcMesh.add_entities(meshset, &tri, 1);
726  CHKERR postProcMesh.add_entities(meshset, edges);
727  CHKERR postProcMesh.convert_entities(meshset, true, false, false);
728  CHKERR postProcMesh.delete_entities(&meshset, 1);
729  CHKERR postProcMesh.delete_entities(edges);
731  };
732 
733  EntityHandle tri;
734 
735  if (elementsMap.find(numeredEntFiniteElementPtr->getEnt()) !=
736  elementsMap.end()) {
737  tri = elementsMap[numeredEntFiniteElementPtr->getEnt()];
738  switch (numeredEntFiniteElementPtr->getEntType()) {
739  case MBTRI:
740  gaussPts.resize(gaussPtsTri.size1(), gaussPtsTri.size2(), false);
741  noalias(gaussPts) = gaussPtsTri;
742  break;
743  case MBQUAD:
744  gaussPts.resize(gaussPtsQuad.size1(), gaussPtsQuad.size2(), false);
745  noalias(gaussPts) = gaussPtsQuad;
746  break;
747  default:
748  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
749  "Not implemented element type");
750  }
751  } else {
752  switch (numeredEntFiniteElementPtr->getEntType()) {
753  case MBTRI:
754  gaussPts.resize(gaussPtsTri.size1(), gaussPtsTri.size2(), false);
755  noalias(gaussPts) = gaussPtsTri;
756  tri = create_tri(gaussPtsTri);
757  break;
758  case MBQUAD:
759  gaussPts.resize(gaussPtsQuad.size1(), gaussPtsQuad.size2(), false);
760  noalias(gaussPts) = gaussPtsQuad;
761  tri = create_quad(gaussPtsQuad);
762  break;
763  default:
764  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
765  "Not implemented element type");
766  }
767  if (sixNodePostProcTris)
768  CHKERR add_mid_nodes(tri);
769  elementsMap[numeredEntFiniteElementPtr->getEnt()] = tri;
770  }
771 
772  // Set values which map nodes with integration points on the prism
773  const EntityHandle *conn;
774  int num_nodes;
775 
776  CHKERR postProcMesh.get_connectivity(tri, conn, num_nodes, false);
777  mapGaussPts.resize(num_nodes);
778  for (int nn = 0; nn != num_nodes; nn++)
779  mapGaussPts[nn] = conn[nn];
780 
782 }
783 
786  ParallelComm *pcomm_post_proc_mesh =
787  ParallelComm::get_pcomm(&postProcMesh, MYPCOMM_INDEX);
788  if (pcomm_post_proc_mesh != NULL) {
789  delete pcomm_post_proc_mesh;
790  }
792 }
793 
796  ParallelComm *pcomm =
797  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
798  ParallelComm *pcomm_post_proc_mesh =
799  ParallelComm::get_pcomm(&postProcMesh, MYPCOMM_INDEX);
800  if (pcomm_post_proc_mesh == NULL) {
801  pcomm_post_proc_mesh = new ParallelComm(&postProcMesh, mField.get_comm());
802  }
803  Range tris;
804  CHKERR postProcMesh.get_entities_by_type(0, MBTRI, tris, false);
805  int rank = pcomm->rank();
806  Range::iterator pit = tris.begin();
807  for (; pit != tris.end(); pit++) {
808  CHKERR postProcMesh.tag_set_data(pcomm_post_proc_mesh->part_tag(), &*pit, 1,
809  &rank);
810  }
811  CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
813 }
814 
817  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
819 
820  if (type != MBVERTEX)
822 
823  CHKERR loopSideVolumes(feVolName, *sideOpFe);
824 
825  // quit if tag is not needed
826  if (!saveOnTag)
828 
829  const MoFEM::FEDofEntity *dof_ptr = data.getFieldDofs()[0].get();
830  int rank = dof_ptr->getNbOfCoeffs();
831 
832  int tag_length = rank * 3;
833  FieldSpace space = dof_ptr->getSpace();
834 
835  double def_VAL[tag_length];
836  bzero(def_VAL, tag_length * sizeof(double));
837  Tag th;
838  CHKERR postProcMesh.tag_get_handle(tagName.c_str(), tag_length,
839  MB_TYPE_DOUBLE, th,
840  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
841 
842  // zero tags, this for Vertex if H1 and TRI if Hdiv, EDGE for Hcurl
843  // no need for L2
844  const void *tags_ptr[mapGaussPts.size()];
845  int nb_gauss_pts = data.getN().size1();
846  if (mapGaussPts.size() != (unsigned int)nb_gauss_pts ||
847  nb_gauss_pts != gradMatPtr->size2()) {
848  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
849  }
850  switch (space) {
851  case H1:
852  case L2:
853 
854  CHKERR postProcMesh.tag_get_by_ptr(th, &mapGaussPts[0], mapGaussPts.size(),
855  tags_ptr);
856  // FIXME: this is not very efficient
857  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
858  for (int rr = 0; rr != rank; ++rr) {
859  for (int dd = 0; dd != 3; ++dd) {
860  const double *my_ptr2 = static_cast<const double *>(tags_ptr[gg]);
861  double *my_ptr = const_cast<double *>(my_ptr2);
862  my_ptr[rank * rr + dd] = (*gradMatPtr)(rank * rr + dd, gg);
863  }
864  }
865  }
866  break;
867  default:
868  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
869  "field with that space is not implemented");
870  }
871 
873 }
874 
876  const std::string field_name, const std::string vol_fe_name,
877  boost::shared_ptr<MatrixDouble> grad_mat_ptr, bool save_on_tag) {
879 
880  if (!grad_mat_ptr)
881  grad_mat_ptr = boost::make_shared<MatrixDouble>();
882 
883  boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> my_side_fe =
884  boost::make_shared<VolumeElementForcesAndSourcesCoreOnSide>(mField);
885 
886  // check number of coefficients
887  auto field_ptr = mField.get_field_structure(field_name);
888  const int nb_coefficients = field_ptr->getNbOfCoeffs();
889 
890  switch (nb_coefficients) {
891  case 1:
892  my_side_fe->getOpPtrVector().push_back(
893  new OpCalculateScalarFieldGradient<3>(field_name, grad_mat_ptr));
894  break;
895  case 3:
896  my_side_fe->getOpPtrVector().push_back(
897  new OpCalculateVectorFieldGradient<3, 3>(field_name, grad_mat_ptr));
898  break;
899  default:
900  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
901  "field with that number of coefficients is not implemented");
902  }
903 
904  FaceElementForcesAndSourcesCore::getOpPtrVector().push_back(
906  postProcMesh, mapGaussPts, field_name, field_name + "_GRAD",
907  my_side_fe, vol_fe_name, grad_mat_ptr, save_on_tag));
908 
910 }
911 
913  return OpSwitch<
914  FaceElementForcesAndSourcesCore::NO_CONTRAVARIANT_TRANSFORM_HDIV |
915  FaceElementForcesAndSourcesCore::NO_COVARIANT_TRANSFORM_HCURL>();
916 }
917 
920 
921  gaussPts.resize(3, 2, false);
922  gaussPts.clear();
923  gaussPts(0, 0) = 0;
924  gaussPts(1, 0) = 0;
925  gaussPts(0, 1) = 1;
926  gaussPts(1, 1) = 0;
927  mapGaussPts.resize(gaussPts.size2());
928 
929  moab::Core core_ref;
930  moab::Interface &moab_ref = core_ref;
931  const EntityHandle *conn;
932  int num_nodes;
933  EntityHandle edge_conn[2];
934  MatrixDouble coords(2, 3);
935  for (int gg = 0; gg != 2; gg++) {
936  coords(gg, 0) = gaussPts(0, gg);
937  coords(gg, 1) = gaussPts(1, gg);
938  coords(gg, 2) = 0;
939  CHKERR moab_ref.create_vertex(&coords(gg, 0), edge_conn[gg]);
940  }
941 
942  EntityHandle edge;
943  CHKERR moab_ref.create_element(MBEDGE, edge_conn, 2, edge);
945 }
946 
949  ParallelComm *pcomm_post_proc_mesh =
950  ParallelComm::get_pcomm(&postProcMesh, MYPCOMM_INDEX);
951  if (pcomm_post_proc_mesh != NULL) {
952  delete pcomm_post_proc_mesh;
953  }
955 }
956 
959  ParallelComm *pcomm =
960  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
961  ParallelComm *pcomm_post_proc_mesh =
962  ParallelComm::get_pcomm(&postProcMesh, MYPCOMM_INDEX);
963  if (pcomm_post_proc_mesh == NULL) {
964  pcomm_post_proc_mesh = new ParallelComm(&postProcMesh, mField.get_comm());
965  }
966  Range edges;
967  CHKERR postProcMesh.get_entities_by_type(0, MBEDGE, edges, false);
968  int rank = pcomm->rank();
969  auto set_edges_rank = [&](const auto rank, const auto &edges) {
970  std::vector<EntityHandle> ranks(edges.size(), rank);
971  CHKERR postProcMesh.tag_set_data(pcomm_post_proc_mesh->part_tag(), edges,
972  &*ranks.begin());
973  };
974  set_edges_rank(rank, edges);
975 
976  CHKERR pcomm_post_proc_mesh->resolve_shared_ents(0);
978 }
979 
980 
983  if (gaussPts.size1() == 0) {
984  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
985  "post-process mesh not generated");
986  }
987 
988  const EntityHandle *conn;
989  int num_nodes;
990  EntityHandle edge;
991 
992  if (elementsMap.find(numeredEntFiniteElementPtr->getEnt()) !=
993  elementsMap.end()) {
994  edge = elementsMap[numeredEntFiniteElementPtr->getEnt()];
995  } else {
996  ublas::vector<EntityHandle> edge_conn(2);
997  MatrixDouble coords_edge(2, 3);
998  VectorDouble coords(3);
999  CHKERR mField.get_moab().get_connectivity(
1000  numeredEntFiniteElementPtr->getEnt(), conn, num_nodes, true);
1001  CHKERR mField.get_moab().get_coords(conn, num_nodes, &coords_edge(0, 0));
1002  for (int gg = 0; gg != 2; gg++) {
1003  double ksi = gaussPts(0, gg);
1004  // double eta = gaussPts(1, gg);
1005  double n0 = N_MBEDGE0(ksi);
1006  double n1 = N_MBEDGE1(ksi);
1007 
1008  double x = n0 * coords_edge(0, 0) + n1 * coords_edge(1, 0);
1009  double y = n0 * coords_edge(0, 1) + n1 * coords_edge(1, 1);
1010  double z = n0 * coords_edge(0, 2) + n1 * coords_edge(1, 2);
1011 
1012  coords[0] = x;
1013  coords[1] = y;
1014  coords[2] = z;
1015 
1016  CHKERR postProcMesh.create_vertex(&coords[0], edge_conn[gg]);
1017  }
1018  CHKERR postProcMesh.create_element(MBEDGE, &edge_conn[0], 2, edge);
1019  elementsMap[numeredEntFiniteElementPtr->getEnt()] = edge;
1020  }
1021 
1022  {
1023  CHKERR postProcMesh.get_connectivity(edge, conn, num_nodes, false);
1024  mapGaussPts.resize(num_nodes);
1025  for (int nn = 0; nn < num_nodes; nn++) {
1026  mapGaussPts[nn] = conn[nn];
1027  }
1028  }
1029 
1031 }
#define N_MBQUAD2(x, y)
quad shape function
Definition: fem_tools.h:67
MoFEMErrorCode generateReferenceElementMesh()
MoFEMErrorCode setGaussPtsTrianglesOnly(int order_triangles_only)
field with continuous normal traction
Definition: definitions.h:179
Post-process fields on refined mesh.
MoFEMErrorCode setGaussPtsThroughThickness(int order_thickness)
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:507
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:75
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
#define N_MBEDGE1(x)
edge shape function
Definition: fem_tools.h:80
MoFEMErrorCode postProcess()
function is run at the end of loop
FieldCoefficientsNumber getNbOfCoeffs() const
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:514
MoFEMErrorCode operator()()
MoFEMErrorCode postProcess()
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
keeps information about indexed dofs for the finite element
double cblas_ddot(const int N, const double *X, const int incX, const double *Y, const int incY)
Definition: cblas_ddot.c:12
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:84
MoFEMErrorCode generateReferenceElementMesh()
MoFEMErrorCode generateReferenceElementMesh()
#define N_MBQUAD3(x, y)
quad shape function
Definition: fem_tools.h:68
FieldSpace getSpace() const
#define N_MBEDGE0(x)
edge shape function
Definition: fem_tools.h:79
#define N_MBQUAD0(x, y)
quad shape function
Definition: fem_tools.h:65
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
field with continuous tangents
Definition: definitions.h:178
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition: Types.hpp:101
FieldSpace
approximation spaces
Definition: definitions.h:174
#define CHKERR
Inline error check.
Definition: definitions.h:602
MoFEMErrorCode preProcess()
constexpr int order
MoFEMErrorCode postProcess()
function is run at the end of loop
#define N_MBTRI2(x, y)
triangle shape function
Definition: fem_tools.h:56
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:91
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:291
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
MoFEMErrorCode preProcess()
function is run at the beginning of loop
MoFEMErrorCode preProcess()
function is run at the beginning of loop
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
#define N_MBQUAD1(x, y)
quad shape function
Definition: fem_tools.h:66
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1788
Data on single entity (This is passed as argument to DataOperator::doWork)
MoFEMErrorCode setGaussPts(int order)
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413
continuous field
Definition: definitions.h:177
static const double eps
MoFEMErrorCode addFieldValuesGradientPostProcOnSkin(const std::string field_name, const std::string vol_fe_name, boost::shared_ptr< MatrixDouble > grad_mat_ptr=nullptr, bool save_on_tag=true)
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:73
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
get dofs values
#define N_MBTRI1(x, y)
triangle shape function
Definition: fem_tools.h:55
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
MoFEMErrorCode setGaussPts(int order)
#define N_MBTRI0(x, y)
triangle shape function
Definition: fem_tools.h:54
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
FTensor::Index< 'i', 2 > i
[Common data]
Definition: ContactOps.hpp:26
field with C-1 continuity
Definition: definitions.h:180