v0.14.0
Tools.cpp
Go to the documentation of this file.
1 /** \file Tools.cpp
2  * \brief Auxiliary tools
3  */
4 
5 #include <phg-quadrule/quad.h>
6 
7 namespace MoFEM {
8 
9 MoFEMErrorCode Tools::query_interface(boost::typeindex::type_index type_index,
10  UnknownInterface **iface) const {
11  *iface = const_cast<Tools *>(this);
12  return 0;
13 }
14 
15 double Tools::volumeLengthQuality(const double *coords) {
16  double lrms = 0;
17  for (int dd = 0; dd != 3; dd++) {
18  lrms += pow(coords[0 * 3 + dd] - coords[1 * 3 + dd], 2) +
19  pow(coords[0 * 3 + dd] - coords[2 * 3 + dd], 2) +
20  pow(coords[0 * 3 + dd] - coords[3 * 3 + dd], 2) +
21  pow(coords[1 * 3 + dd] - coords[2 * 3 + dd], 2) +
22  pow(coords[1 * 3 + dd] - coords[3 * 3 + dd], 2) +
23  pow(coords[2 * 3 + dd] - coords[3 * 3 + dd], 2);
24  }
25  lrms = sqrt((1. / 6.) * lrms);
26  double volume = tetVolume(coords);
27  return 6. * sqrt(2.) * volume / pow(lrms, 3);
28 }
29 
30 double Tools::tetVolume(const double *coords) {
31  double diff_n[12];
32  ShapeDiffMBTET(diff_n);
33  FTensor::Tensor1<double *, 3> t_diff_n(&diff_n[0], &diff_n[1], &diff_n[2], 3);
34  FTensor::Tensor1<const double *, 3> t_coords(&coords[0], &coords[1],
35  &coords[2], 3);
36  FTensor::Tensor2<double, 3, 3> jac;
37  FTensor::Index<'i', 3> i;
38  FTensor::Index<'j', 3> j;
39  jac(i, j) = 0;
40  for (int nn = 0; nn != 4; nn++) {
41  jac(i, j) += t_coords(i) * t_diff_n(j);
42  ++t_coords;
43  ++t_diff_n;
44  }
45  return determinantTensor3by3(jac) / 6.;
46 }
47 
48 MoFEMErrorCode
49 Tools::minTetsQuality(const Range &tets, double &min_quality, Tag th,
50  boost::function<double(double, double)> f) {
51  MoFEM::Interface &m_field = cOre;
52  moab::Interface &moab(m_field.get_moab());
53  MoFEMFunctionBegin;
54  const EntityHandle *conn;
55  int num_nodes;
56  double coords[12];
57  for (auto tet : tets) {
58  CHKERR m_field.get_moab().get_connectivity(tet, conn, num_nodes, true);
59  if (th) {
60  CHKERR moab.tag_get_data(th, conn, num_nodes, coords);
61  } else {
62  CHKERR moab.get_coords(conn, num_nodes, coords);
63  }
64  double q = Tools::volumeLengthQuality(coords);
65  if (!std::isnormal(q))
66  q = -2;
67  min_quality = f(q, min_quality);
68  }
69  MoFEMFunctionReturn(0);
70 }
71 
72 constexpr double Tools::shapeFunMBEDGE0At00;
73 constexpr double Tools::shapeFunMBEDGE1At00;
74 constexpr std::array<double, 2> Tools::shapeFunMBEDGEAt00;
75 
76 constexpr std::array<double, 2> Tools::diffShapeFunMBEDGE;
77 constexpr std::array<double, 6> Tools::diffShapeFunMBTRI;
78 constexpr std::array<double, 12> Tools::diffShapeFunMBTET;
79 constexpr double Tools::shapeFunMBTRI0At00;
80 constexpr double Tools::shapeFunMBTRI1At00;
81 constexpr double Tools::shapeFunMBTRI2At00;
82 constexpr std::array<double, 3> Tools::shapeFunMBTRIAt00;
83 
84 constexpr std::array<double, 4> Tools::shapeFunMBTETAt000;
85 constexpr std::array<double, 8> Tools::diffShapeFunMBQUADAtCenter;
86 constexpr std::array<double, 24> Tools::diffShapeFunMBHEXAtCenter;
87 
88 MoFEMErrorCode Tools::getLocalCoordinatesOnReferenceFourNodeTet(
89  const double *elem_coords, const double *global_coords, const int nb_nodes,
90  double *local_coords) {
91  FTensor::Index<'i', 4> i;
92  MoFEMFunctionBeginHot;
93 
94  FTensor::Tensor1<FTensor::PackPtr<const double *, 1>, 4> t_elem_coords = {
95  &elem_coords[0], &elem_coords[3], &elem_coords[6], &elem_coords[9]};
96 
97  FTensor::Tensor1<const double, 4> t_n = {
98  shapeFunMBTETAt000[0], shapeFunMBTETAt000[1], shapeFunMBTETAt000[2],
99  shapeFunMBTETAt000[3]};
100  FTensor::Tensor1<double, 3> t_coords_at_0;
101 
102  // Build matrix and get coordinates of zero point
103  // ii - global coordinates
104  // jj - local derivatives
105  MatrixDouble3by3 a(3, 3);
106  for (auto ii : {0, 1, 2}) {
107  FTensor::Tensor1<FTensor::PackPtr<const double *, 1>, 4> t_diff(
108  &diffShapeFunMBTET[0], &diffShapeFunMBTET[3], &diffShapeFunMBTET[6],
109  &diffShapeFunMBTET[9]);
110  for (auto jj : {0, 1, 2}) {
111  a(jj, ii) = t_diff(i) * t_elem_coords(i);
112  ++t_diff;
113  }
114  t_coords_at_0(ii) = t_n(i) * t_elem_coords(i);
115  ++t_elem_coords;
116  }
117 
118  FTensor::Tensor1<FTensor::PackPtr<const double *, 3>, 3> t_global_coords = {
119  &global_coords[0], &global_coords[1], &global_coords[2]};
120  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_local_coords = {
121  &local_coords[0], &local_coords[1], &local_coords[2]};
122 
123  // Calculate right hand side
124  FTensor::Index<'j', 3> j;
125  for (int ii = 0; ii != nb_nodes; ++ii) {
126  t_local_coords(j) = t_global_coords(j) - t_coords_at_0(j);
127  ++t_local_coords;
128  ++t_global_coords;
129  }
130 
131  // Solve problem
132  int IPIV[3];
133  int info = lapack_dgesv(3, nb_nodes, &a(0, 0), 3, IPIV, local_coords, 3);
134  if (info != 0)
135  SETERRQ1(PETSC_COMM_SELF, 1, "info == %d", info);
136 
137  MoFEMFunctionReturnHot(0);
138 }
139 
140 template <typename T1, typename T2>
141 MoFEMErrorCode getLocalCoordinatesOnReferenceThreeNodeTriImpl(
142  const T1 *elem_coords, const T2 *global_coords, const int nb_nodes,
143  typename FTensor::promote<T1, T2>::V *local_coords) {
144 
145  FTensor::Index<'i', 3> i3;
146  FTensor::Index<'j', 3> j3;
147  FTensor::Index<'K', 2> K2;
148  FTensor::Index<'L', 2> L2;
149  MoFEMFunctionBeginHot;
150 
151  FTensor::Tensor1<const double, 3> t_n = {Tools::shapeFunMBTRIAt00[0],
152  Tools::shapeFunMBTRIAt00[1],
153  Tools::shapeFunMBTRIAt00[2]};
154  auto t_diff = getFTensor2FromPtr<3, 2>(
155  const_cast<double *>(Tools::diffShapeFunMBTRI.data()));
156  auto t_elem_coords = getFTensor2FromPtr<3, 3>(const_cast<T1 *>(elem_coords));
157 
158  // Build matrix and get coordinates of zero point
159  FTensor::Tensor2<T1, 2, 3> t_a;
160  t_a(K2, i3) = t_diff(j3, K2) * t_elem_coords(j3, i3);
161  FTensor::Tensor2_symmetric<T1, 2> t_b, t_inv_b;
162  t_b(K2, L2) = t_a(K2, j3) ^ t_a(L2, j3);
163  // Solve problem
164  const auto inv_det = 1. / (t_b(0, 0) * t_b(1, 1) - t_b(0, 1) * t_b(1, 0));
165  t_inv_b(0, 0) = t_b(1, 1) * inv_det;
166  t_inv_b(0, 1) = -t_b(0, 1) * inv_det;
167  t_inv_b(1, 1) = t_b(0, 0) * inv_det;
168 
169  // Coords at corner
170  FTensor::Tensor1<T1, 3> t_coords_at_0;
171  t_coords_at_0(i3) = t_n(j3) * t_elem_coords(j3, i3);
172 
173  auto t_global_coords = getFTensor1FromPtr<3>(const_cast<T2 *>(global_coords));
174  auto t_local_coords = getFTensor1FromPtr<2>(local_coords);
175 
176  // Calculate right hand side
177  for (int ii = 0; ii != nb_nodes; ++ii) {
178  t_local_coords(L2) =
179  t_inv_b(L2, K2) *
180  (t_a(K2, j3) * (t_global_coords(j3) - t_coords_at_0(j3)));
181  ++t_local_coords;
182  ++t_global_coords;
183  }
184 
185  MoFEMFunctionReturnHot(0);
186 }
187 
188 MoFEMErrorCode Tools::getLocalCoordinatesOnReferenceThreeNodeTri(
189  const double *elem_coords, const double *glob_coords, const int nb_nodes,
190  double *local_coords) {
191  return getLocalCoordinatesOnReferenceThreeNodeTriImpl<double, double>(
192  elem_coords, glob_coords, nb_nodes, local_coords);
193 }
194 
195 MoFEMErrorCode Tools::getLocalCoordinatesOnReferenceThreeNodeTri(
196  const double *elem_coords, const std::complex<double> *glob_coords,
197  const int nb_nodes, std::complex<double> *local_coords) {
198  return getLocalCoordinatesOnReferenceThreeNodeTriImpl<double,
199  std::complex<double>>(
200  elem_coords, glob_coords, nb_nodes, local_coords);
201 }
202 
203 MoFEMErrorCode Tools::getLocalCoordinatesOnReferenceEdgeNodeEdge(
204  const double *elem_coords, const double *global_coords, const int nb_nodes,
205  double *local_coords) {
206 
207  FTensor::Index<'i', 3> i3;
208  MoFEMFunctionBeginHot;
209 
210  FTensor::Tensor1<FTensor::PackPtr<const double *, 1>, 3> t_elem_coords = {
211  &elem_coords[0], &elem_coords[3], &elem_coords[6]};
212 
213  FTensor::Tensor1<double, 3> t_coords_at_0;
214  // Build matrix and get coordinates of zero point
215  // ii - global coordinates
216  FTensor::Tensor1<double, 3> t_a;
217  for (auto ii : {0, 1, 2}) {
218  t_a(ii) = diffShapeFunMBEDGE[0] * t_elem_coords(0) +
219  diffShapeFunMBEDGE[1] * t_elem_coords(1);
220  t_coords_at_0(ii) = shapeFunMBEDGEAt00[0] * t_elem_coords(0) +
221  shapeFunMBEDGEAt00[1] * t_elem_coords(1);
222  ++t_elem_coords;
223  }
224 
225  FTensor::Tensor1<FTensor::PackPtr<const double *, 3>, 3> t_global_coords = {
226  &global_coords[0], &global_coords[1], &global_coords[2]};
227  FTensor::Tensor0<FTensor::PackPtr<double *, 1>> t_local_coords = {
228  &local_coords[0]};
229 
230  const double b = t_a(i3) * t_a(i3);
231  const double inv_b = 1 / b;
232 
233  // Calculate right hand side
234  for (int ii = 0; ii != nb_nodes; ++ii) {
235  t_local_coords =
236  inv_b * (t_a(i3) * (t_global_coords(i3) - t_coords_at_0(i3)));
237  ++t_local_coords;
238  ++t_global_coords;
239  }
240 
241  MoFEMFunctionReturnHot(0);
242 }
243 
244 MoFEMErrorCode Tools::getTetsWithQuality(Range &out_tets, const Range &tets,
245  Tag th,
246  boost::function<bool(double)> f) {
247  MoFEM::Interface &m_field = cOre;
248  moab::Interface &moab(m_field.get_moab());
249  MoFEMFunctionBegin;
250  Range to_write;
251  const EntityHandle *conn;
252  int num_nodes;
253  double coords[12];
254  for (auto tet : tets) {
255  CHKERR m_field.get_moab().get_connectivity(tet, conn, num_nodes, true);
256  if (th) {
257  CHKERR moab.tag_get_data(th, conn, num_nodes, coords);
258  } else {
259  CHKERR moab.get_coords(conn, num_nodes, coords);
260  }
261  double q = Tools::volumeLengthQuality(coords);
262  if (f(q)) {
263  out_tets.insert(tet);
264  }
265  }
266  MoFEMFunctionReturn(0);
267 }
268 
269 MoFEMErrorCode Tools::writeTetsWithQuality(const char *file_name,
270  const char *file_type,
271  const char *options,
272  const Range &tets, Tag th,
273  boost::function<bool(double)> f) {
274  MoFEM::Interface &m_field = cOre;
275  moab::Interface &moab(m_field.get_moab());
276  MoFEMFunctionBegin;
277  Range out_tets;
278  CHKERR getTetsWithQuality(out_tets, tets, th, f);
279  EntityHandle meshset;
280  CHKERR moab.create_meshset(MESHSET_SET, meshset);
281  CHKERR moab.add_entities(meshset, out_tets);
282  CHKERR moab.write_file(file_name, file_type, options, &meshset, 1);
283  CHKERR moab.delete_entities(&meshset, 1);
284  MoFEMFunctionReturn(0);
285 }
286 
287 MoFEMErrorCode Tools::checkIfPointIsInTet(const double tet_coords[],
288  const double global_coord[],
289  const double tol, bool &result) {
290  double loc_coord[] = {0, 0, 0};
291  double N[4], diffN[12];
292  MoFEMFunctionBegin;
293  CHKERR ShapeDiffMBTET(diffN);
294  CHKERR ShapeMBTET(N, &loc_coord[0], &loc_coord[1], &loc_coord[2], 1);
295  CHKERR ShapeMBTET_inverse(N, diffN, tet_coords, global_coord, loc_coord);
296  CHKERR ShapeMBTET(N, &loc_coord[0], &loc_coord[1], &loc_coord[2], 1);
297  result = true;
298  for (int n = 0; n != 4; ++n) {
299  if (N[n] < -tol || (N[n] - 1) > tol) {
300  result = false;
301  break;
302  }
303  }
304  MoFEMFunctionReturn(0);
305 }
306 
307 MoFEMErrorCode Tools::checkVectorForNotANumber(const Problem *prb_ptr,
308  const RowColData row_or_col,
309  Vec v) {
310  MoFEMFunctionBegin;
311  int loc_size;
312  CHKERR VecGetLocalSize(v, &loc_size);
313  int prb_loc_size = 0;
314  boost::shared_ptr<NumeredDofEntity_multiIndex> prb_dofs;
315  switch (row_or_col) {
316  case ROW:
317  prb_loc_size = prb_ptr->getNbLocalDofsRow();
318  prb_dofs = prb_ptr->getNumeredRowDofsPtr();
319  break;
320  case COL:
321  prb_loc_size = prb_ptr->getNbLocalDofsCol();
322  prb_dofs = prb_ptr->getNumeredColDofsPtr();
323  break;
324  break;
325  default:
326  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
327  "Wrong argument, row_or_col should be row or column");
328  }
329  if (loc_size != prb_loc_size) {
330  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
331  "Inconsistent size of vector and problem %d != %d", loc_size,
332  prb_loc_size);
333  }
334  const double *a;
335  CHKERR VecGetArrayRead(v, &a);
336  MPI_Comm comm = PetscObjectComm((PetscObject)v);
337  for (int ii = 0; ii != loc_size; ++ii) {
338  if (!boost::math::isfinite(a[ii])) {
339  NumeredDofEntityByLocalIdx::iterator dit =
340  prb_dofs->get<PetscLocalIdx_mi_tag>().find(ii);
341  std::ostringstream ss;
342  ss << "Not a number " << a[ii] << " on dof: " << endl
343  << **dit << endl
344  << endl;
345  PetscSynchronizedPrintf(comm, "%s", ss.str().c_str());
346  }
347  }
348  CHKERR VecRestoreArrayRead(v, &a);
349  PetscSynchronizedFlush(comm, PETSC_STDOUT);
350  MoFEMFunctionReturn(0);
351 }
352 
353 MoFEMErrorCode Tools::getTriNormal(const double *coords, double *normal,
354  double *d_normal) {
355  MoFEMFunctionBegin;
356  FTensor::Index<'i', 3> i;
357  FTensor::Index<'j', 3> j;
358  FTensor::Index<'k', 3> k;
359  FTensor::Index<'l', 3> l;
360  FTensor::Index<'n', 3> n;
361  FTensor::Index<'m', 3> m;
362  FTensor::Index<'J', 2> J;
363  FTensor::Number<0> N0;
364  FTensor::Number<1> N1;
365  auto diff_ptr = Tools::diffShapeFunMBTRI.data();
366  auto t_diff_tensor = getFTensor2FromPtr<3, 2>(const_cast<double *>(diff_ptr));
367  auto t_coords = getFTensor2FromPtr<3, 3>(const_cast<double *>(coords));
368  FTensor::Tensor2<double, 3, 2> t_tangent;
369  t_tangent(i, J) = t_coords(n, i) * t_diff_tensor(n, J);
370  auto t_normal = getFTensor1FromPtr<3>(normal);
371  t_normal(j) =
372  FTensor::levi_civita(i, j, k) * t_tangent(k, N0) * t_tangent(i, N1);
373  if (d_normal) {
374  constexpr auto t_kd = FTensor::Kronecker_Delta<int>();
375  FTensor::Tensor4<int, 3, 3, 3, 3> t_d_coords;
376  t_d_coords(i, j, k, n) = t_kd(i, k) * t_kd(j, n);
377  FTensor::Tensor4<double, 3, 3, 3, 2> t_d_tangent;
378  t_d_tangent(i, k, l, J) = t_d_coords(n, i, k, l) * t_diff_tensor(n, J);
379  auto t_d_normal = getFTensor3FromPtr<3, 3, 3>(d_normal);
380  t_d_normal(j, m, n) = (FTensor::levi_civita(i, j, k) * t_tangent(i, N1)) *
381  t_d_tangent(k, m, n, N0)
382 
383  +
384 
385  (FTensor::levi_civita(i, j, k) * t_tangent(k, N0)) *
386  t_d_tangent(i, m, n, N1);
387  }
388  MoFEMFunctionReturn(0);
389 }
390 
391 MoFEMErrorCode Tools::getTriNormal(const EntityHandle tri,
392  double *normal) const {
393  MoFEM::Interface &m_field = cOre;
394  moab::Interface &moab(m_field.get_moab());
395  MoFEMFunctionBegin;
396  if (type_from_handle(tri) != MBTRI) {
397  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "Works only for triangle");
398  }
399  const EntityHandle *conn;
400  int num_nodes;
401  double coords[9];
402  CHKERR moab.get_connectivity(tri, conn, num_nodes, true);
403  CHKERR moab.get_coords(conn, num_nodes, coords);
404  CHKERR getTriNormal(coords, normal);
405  MoFEMFunctionReturn(0);
406 }
407 
408 double Tools::getTriArea(const EntityHandle tri) const {
409  FTensor::Index<'i', 3> i;
410  FTensor::Tensor1<double, 3> t_normal;
411  CHK_THROW_MESSAGE(getTriNormal(tri, &t_normal(0)), "calculate area");
412  return sqrt(t_normal(i) * t_normal(i)) * 0.5;
413 }
414 
415 double Tools::getEdgeLength(const double *edge_coords) {
416  FTensor::Tensor1<double, 3> t_coords_n0(edge_coords[0], edge_coords[1],
417  edge_coords[2]);
418  FTensor::Tensor1<double, 3> t_coords_n1(edge_coords[3], edge_coords[4],
419  edge_coords[5]);
420  FTensor::Index<'i', 3> i;
421  t_coords_n0(i) -= t_coords_n1(i);
422  return sqrt(t_coords_n0(i) * t_coords_n0(i));
423 }
424 
425 double Tools::getEdgeLength(const EntityHandle edge) {
426  MoFEM::Interface &m_field = cOre;
427  moab::Interface &moab(m_field.get_moab());
428  auto get_edge_coords = [edge, &moab](double *const coords) {
429  MoFEMFunctionBegin;
430  if (type_from_handle(edge) != MBEDGE) {
431  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "Works only for edge");
432  }
433  const EntityHandle *conn;
434  int num_nodes;
435  CHKERR moab.get_connectivity(edge, conn, num_nodes, true);
436  CHKERR moab.get_coords(conn, 2, coords);
437  MoFEMFunctionReturn(0);
438  };
439  double coords[6];
440  ierr = get_edge_coords(coords);
441  CHKERRABORT(PETSC_COMM_SELF, ierr);
442  return getEdgeLength(coords);
443 }
444 
445 Tools::SEGMENT_MIN_DISTANCE
446 Tools::minDistancePointFromOnSegment(const double *w_ptr, const double *v_ptr,
447  const double *p_ptr, double *const t_ptr) {
448  FTensor::Index<'i', 3> i;
449  FTensor::Tensor1<const double *, 3> t_w(&w_ptr[0], &w_ptr[1], &w_ptr[2]);
450  FTensor::Tensor1<const double *, 3> t_v(&v_ptr[0], &v_ptr[1], &v_ptr[2]);
451  FTensor::Tensor1<const double *, 3> t_p(&p_ptr[0], &p_ptr[1], &p_ptr[2]);
452  FTensor::Tensor1<double, 3> t_vw;
453  t_vw(i) = t_v(i) - t_w(i);
454  const double dot_vw = t_vw(i) * t_vw(i);
455  if (std::abs(dot_vw) < std::numeric_limits<double>::epsilon()) {
456  if (t_ptr)
457  *t_ptr = 0;
458  return SEGMENT_ONE_IS_POINT;
459  }
460  FTensor::Tensor1<double, 3> t_pw;
461  t_pw(i) = t_p(i) - t_w(i);
462  const double t = t_pw(i) * t_vw(i) / dot_vw;
463  if (t_ptr)
464  *t_ptr = t;
465  return SOLUTION_EXIST;
466 }
467 
468 Tools::SEGMENT_MIN_DISTANCE
469 Tools::minDistanceFromSegments(const double *w_ptr, const double *v_ptr,
470  const double *k_ptr, const double *l_ptr,
471  double *const tvw_ptr, double *const tlk_ptr) {
472 
473  FTensor::Index<'i', 3> i;
474  FTensor::Tensor1<const double *, 3> t_w(&w_ptr[0], &w_ptr[1], &w_ptr[2]);
475  FTensor::Tensor1<const double *, 3> t_v(&v_ptr[0], &v_ptr[1], &v_ptr[2]);
476  FTensor::Tensor1<const double *, 3> t_k(&k_ptr[0], &k_ptr[1], &k_ptr[2]);
477  FTensor::Tensor1<const double *, 3> t_l(&l_ptr[0], &l_ptr[1], &l_ptr[2]);
478 
479  // First segnent is a point
480  FTensor::Tensor1<double, 3> t_vw;
481  t_vw(i) = t_v(i) - t_w(i);
482  double dot_vw = t_vw(i) * t_vw(i);
483  if (std::abs(dot_vw) < std::numeric_limits<double>::epsilon()) {
484  if (tvw_ptr)
485  *tvw_ptr = 0;
486  if (minDistancePointFromOnSegment(k_ptr, l_ptr, w_ptr, tlk_ptr) ==
487  SEGMENT_ONE_IS_POINT)
488  return SEGMENT_TWO_AND_TWO_ARE_POINT;
489  else
490  return SEGMENT_ONE_IS_POINT;
491  }
492 
493  // Second segment is a point
494  FTensor::Tensor1<double, 3> t_lk;
495  t_lk(i) = t_l(i) - t_k(i);
496  double dot_lk = t_lk(i) * t_lk(i);
497  if (std::abs(dot_lk) < std::numeric_limits<double>::epsilon()) {
498  if (tlk_ptr)
499  *tlk_ptr = 0;
500  if (minDistancePointFromOnSegment(w_ptr, v_ptr, k_ptr, tvw_ptr) ==
501  SEGMENT_ONE_IS_POINT)
502  return SEGMENT_TWO_AND_TWO_ARE_POINT;
503  else
504  return SEGMENT_TWO_IS_POINT;
505  }
506 
507  const double a = t_vw(i) * t_vw(i);
508  const double b = -t_vw(i) * t_lk(i);
509  const double c = t_lk(i) * t_lk(i);
510 
511  const double det = a * c - b * b;
512  if (std::abs(det) < std::numeric_limits<double>::epsilon()) {
513 
514  return NO_SOLUTION;
515 
516  } else {
517 
518  FTensor::Tensor1<double, 3> t_wk;
519  t_wk(i) = t_w(i) - t_k(i);
520 
521  const double ft0 = t_vw(i) * t_wk(i);
522  const double ft1 = -t_lk(i) * t_wk(i);
523  const double t0 = (ft1 * b - ft0 * c) / det;
524  const double t1 = (ft0 * b - ft1 * a) / det;
525 
526  if (tvw_ptr)
527  *tvw_ptr = t0;
528  if (tlk_ptr)
529  *tlk_ptr = t1;
530 
531  return SOLUTION_EXIST;
532  }
533 }
534 
535 MoFEMErrorCode Tools::findMinDistanceFromTheEdges(
536  const double *v_ptr, const int nb, Range edges, double *min_dist_ptr,
537  double *o_ptr, EntityHandle *o_segments) const {
538  MoFEM::Interface &m_field = cOre;
539  moab::Interface &moab(m_field.get_moab());
540  MoFEMFunctionBegin;
541 
542  FTensor::Index<'i', 3> i;
543 
544  auto get_point = [i](auto &t_w, auto &t_delta, auto t) {
545  FTensor::Tensor1<double, 3> t_p;
546  t = std::max(0., std::min(1., t));
547  t_p(i) = t_w(i) + t * t_delta(i);
548  return t_p;
549  };
550 
551  auto get_distance = [i](auto &t_p, auto &t_n) {
552  FTensor::Tensor1<double, 3> t_dist_vector;
553  t_dist_vector(i) = t_p(i) - t_n(i);
554  return sqrt(t_dist_vector(i) * t_dist_vector(i));
555  };
556 
557  for (auto e : edges) {
558 
559  int num_nodes;
560  const EntityHandle *conn_fixed;
561  CHKERR moab.get_connectivity(e, conn_fixed, num_nodes, true);
562  VectorDouble6 coords_fixed(6);
563  CHKERR moab.get_coords(conn_fixed, num_nodes, &coords_fixed[0]);
564  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_f0(
565  &coords_fixed[0], &coords_fixed[1], &coords_fixed[2]);
566  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_f1(
567  &coords_fixed[3], &coords_fixed[4], &coords_fixed[5]);
568 
569  FTensor::Tensor1<double, 3> t_edge_delta;
570  t_edge_delta(i) = t_f1(i) - t_f0(i);
571 
572  FTensor::Tensor1<FTensor::PackPtr<const double *, 3>, 3> t_n(
573  v_ptr, v_ptr + 1, v_ptr + 2);
574  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_min_coords(
575  o_ptr, o_ptr + 1, o_ptr + 2);
576  FTensor::Tensor0<FTensor::PackPtr<double *, 1>> t_min_dist(min_dist_ptr);
577 
578  EntityHandle *colsest_segment_it = nullptr;
579  if (o_segments)
580  colsest_segment_it = o_segments;
581 
582  for (int n = 0; n != nb; ++n) {
583 
584  double t;
585  if (Tools::minDistancePointFromOnSegment(&t_f0(0), &t_f1(0), &t_n(0),
586  &t) == Tools::SOLUTION_EXIST) {
587  auto t_p = get_point(t_f0, t_edge_delta, t);
588  auto dist_n = get_distance(t_p, t_n);
589  if (dist_n < t_min_dist || t_min_dist < 0) {
590  t_min_dist = dist_n;
591  if (o_ptr)
592  t_min_coords(i) = t_p(i);
593  if (o_segments)
594  *colsest_segment_it = e;
595  }
596  }
597 
598  ++t_n;
599  ++t_min_dist;
600  if (o_ptr)
601  ++t_min_coords;
602  if (o_segments)
603  ++colsest_segment_it;
604  }
605  }
606 
607  MoFEMFunctionReturn(0);
608 }
609 
610 MoFEMErrorCode Tools::outerProductOfEdgeIntegrationPtsForQuad(
611  MatrixDouble &gauss_pts, const int rule_ksi, const int rule_eta) {
612  MoFEMFunctionBegin;
613 
614  auto check_rule_edge = [](int rule) {
615  MoFEMFunctionBeginHot;
616  if (rule < 0) {
617  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
618  "Wrong integration rule: %d", rule);
619  }
620  if (rule > QUAD_1D_TABLE_SIZE) {
621  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
622  "rule > quadrature order %d < %d", rule, QUAD_1D_TABLE_SIZE);
623  }
624  if (QUAD_1D_TABLE[rule]->dim != 1) {
625  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong dimension");
626  }
627  if (QUAD_1D_TABLE[rule]->order < rule) {
628  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
629  "wrong order %d != %d", QUAD_1D_TABLE[rule]->order, rule);
630  }
631  MoFEMFunctionReturnHot(0);
632  };
633 
634  CHKERR check_rule_edge(rule_ksi);
635  CHKERR check_rule_edge(rule_eta);
636 
637  const int nb_gauss_pts_ksi = QUAD_1D_TABLE[rule_ksi]->npoints;
638  const int nb_gauss_pts_eta = QUAD_1D_TABLE[rule_eta]->npoints;
639  gauss_pts.resize(3, nb_gauss_pts_ksi * nb_gauss_pts_eta, false);
640 
641  int gg = 0;
642  for (size_t i = 0; i != nb_gauss_pts_ksi; ++i) {
643  const double wi = QUAD_1D_TABLE[rule_ksi]->weights[i];
644  const double ksi = (QUAD_1D_TABLE[rule_ksi]->points[2 * i + 1]);
645  for (size_t j = 0; j != nb_gauss_pts_eta; ++j, ++gg) {
646  const double wk = wi * QUAD_1D_TABLE[rule_eta]->weights[j];
647  const double eta = QUAD_1D_TABLE[rule_eta]->points[2 * j + 1];
648  gauss_pts(0, gg) = ksi;
649  gauss_pts(1, gg) = eta;
650  gauss_pts(2, gg) = wk;
651  }
652  }
653  if (gg != gauss_pts.size2())
654  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong size of matrix");
655 
656  MoFEMFunctionReturn(0);
657 }
658 
659 MoFEMErrorCode Tools::outerProductOfEdgeIntegrationPtsForHex(
660  MatrixDouble &gauss_pts, const int rule_ksi, const int rule_eta,
661  const int rule_zeta) {
662  MoFEMFunctionBegin;
663 
664  auto check_rule_edge = [](int rule) {
665  MoFEMFunctionBeginHot;
666  if (rule < 0) {
667  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
668  "Wrong integration rule: %d", rule);
669  }
670  if (rule > QUAD_1D_TABLE_SIZE) {
671  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
672  "rule > quadrature order %d < %d", rule, QUAD_1D_TABLE_SIZE);
673  }
674  if (QUAD_1D_TABLE[rule]->dim != 1) {
675  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong dimension");
676  }
677  if (QUAD_1D_TABLE[rule]->order < rule) {
678  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
679  "wrong order %d != %d", QUAD_1D_TABLE[rule]->order, rule);
680  }
681  MoFEMFunctionReturnHot(0);
682  };
683 
684  CHKERR check_rule_edge(rule_ksi);
685  CHKERR check_rule_edge(rule_eta);
686  CHKERR check_rule_edge(rule_zeta);
687 
688  const int nb_gauss_pts_ksi = QUAD_1D_TABLE[rule_ksi]->npoints;
689  const int nb_gauss_pts_eta = QUAD_1D_TABLE[rule_eta]->npoints;
690  const int nb_gauss_pts_zeta = QUAD_1D_TABLE[rule_zeta]->npoints;
691  gauss_pts.resize(4, nb_gauss_pts_ksi * nb_gauss_pts_eta * nb_gauss_pts_zeta,
692  false);
693 
694  int gg = 0;
695  for (size_t i = 0; i != nb_gauss_pts_ksi; ++i) {
696  const double wi = QUAD_1D_TABLE[rule_ksi]->weights[i];
697  const double ksi = QUAD_1D_TABLE[rule_ksi]->points[2 * i + 1];
698  for (size_t j = 0; j != nb_gauss_pts_eta; ++j) {
699  const double wj = wi * QUAD_1D_TABLE[rule_eta]->weights[j];
700  const double eta = QUAD_1D_TABLE[rule_eta]->points[2 * j + 1];
701  for (size_t k = 0; k != nb_gauss_pts_zeta; ++k, ++gg) {
702  const double wk = wj * QUAD_1D_TABLE[rule_zeta]->weights[k];
703  const double zeta = QUAD_1D_TABLE[rule_zeta]->points[2 * k + 1];
704  gauss_pts(0, gg) = ksi;
705  gauss_pts(1, gg) = eta;
706  gauss_pts(2, gg) = zeta;
707  gauss_pts(3, gg) = wk;
708  }
709  }
710  }
711 
712  if (gg != gauss_pts.size2())
713  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong size of matrix");
714 
715  MoFEMFunctionReturn(0);
716 }
717 
718 constexpr std::array<int, 12> Tools::uniformTriangleRefineTriangles;
719 
720 std::tuple<std::vector<double>, std::vector<int>, std::vector<int>>
721 Tools::refineTriangle(int nb_levels) {
722 
723  std::vector<int> triangles{0, 1, 2, 3, 4, 5};
724  std::vector<double> nodes{
725 
726  0., 0., // 0
727  1., 0., // 1
728  0., 1., // 2
729  0.5, 0., // 3
730  0.5, 0.5, // 4
731  0., 0.5 // 5
732 
733  };
734  std::map<std::pair<int, int>, int> edges{
735  {{0, 1}, 3}, {{1, 2}, 4}, {{0, 2}, 5}};
736 
737  auto add_edge = [&](auto a, auto b) {
738  if (a > b) {
739  std::swap(a, b);
740  }
741  auto it = edges.find(std::make_pair(a, b));
742  if (it == edges.end()) {
743  int e = 3 + edges.size();
744  edges[std::make_pair(a, b)] = e;
745  for (auto n : {0, 1}) {
746  nodes.push_back((nodes[2 * a + n] + nodes[2 * b + n]) / 2);
747  }
748 #ifndef NDEBUG
749  if (e != nodes.size() / 2 - 1)
750  CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "wrong node/edge index");
751 #endif
752  return e;
753  }
754  return it->second;
755  };
756 
757  auto add_triangle = [&](auto t) {
758  for (auto tt : {0, 1, 2, 3}) {
759  auto last = triangles.size() / 6;
760  for (auto n : {0, 1, 2}) {
761  // add triangle nodes
762  triangles.push_back(
763  triangles[6 * t + uniformTriangleRefineTriangles[3 * tt + n]]);
764  }
765  // add triangle edges
766  auto cycle = std::array<int, 4>{0, 1, 2, 0};
767  for (auto e : {0, 1, 2}) {
768  triangles.push_back(add_edge(triangles[6 * last + cycle[e]],
769  triangles[6 * last + cycle[e + 1]]));
770  }
771  }
772  };
773 
774  std::vector<int> level_index{0, 1};
775  auto l = 0;
776  for (; l != nb_levels; ++l) {
777  auto first_tet = level_index[l];
778  auto nb_last_level_test = level_index.back() - level_index[l];
779  for (auto t = first_tet; t != (first_tet + nb_last_level_test); ++t) {
780  add_triangle(t);
781  }
782  level_index.push_back(triangles.size() / 6);
783  }
784 
785  return std::make_tuple(nodes, triangles, level_index);
786 }
787 
788 MatrixDouble Tools::refineTriangleIntegrationPts(
789  MatrixDouble pts,
790  std::tuple<std::vector<double>, std::vector<int>, std::vector<int>>
791  refined) {
792 
793  auto [nodes, triangles, level_index] = refined;
794 
795  auto get_coords = [&](auto t) {
796  auto [nodes, triangles, level_index] = refined;
797  std::array<double, 9> ele_coords;
798  for (auto n : {0, 1, 2}) {
799  for (auto i : {0, 1}) {
800  ele_coords[3 * n + i] = nodes[2 * triangles[6 * t + n] + i];
801  }
802  ele_coords[3 * n + 2] = 0;
803  }
804  return ele_coords;
805  };
806 
807  auto get_normal = [](auto &ele_coords) {
808  FTensor::Tensor1<double, 3> t_normal;
809  Tools::getTriNormal(ele_coords.data(), &t_normal(0));
810  return t_normal;
811  };
812 
813  std::vector<double> ele_shape(3 * pts.size2());
814  shapeFunMBTRI<1>(&*ele_shape.begin(), &pts(0, 0), &pts(1, 0), pts.size2());
815 
816  int nb_elems = level_index.back() - level_index[level_index.size() - 2];
817  MatrixDouble new_pts(3, pts.size2() * nb_elems);
818  FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 2> t_gauss_pt{&new_pts(0, 0),
819  &new_pts(1, 0)};
820  FTensor::Tensor0<FTensor::PackPtr<double *, 1>> t_gauss_weight{
821  &new_pts(2, 0)};
822 
823  for (auto t = level_index[level_index.size() - 2]; t != level_index.back();
824  ++t) {
825 
826  auto ele_coords = get_coords(t);
827  auto t_normal = get_normal(ele_coords);
828  auto area = t_normal.l2();
829 
830  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_ele_shape{
831  &ele_shape[0], &ele_shape[1], &ele_shape[2]};
832  FTensor::Tensor2<double, 3, 2> t_ele_coords{ele_coords[0], ele_coords[1],
833  ele_coords[3], ele_coords[4],
834  ele_coords[6], ele_coords[7]};
835 
836  FTensor::Index<'i', 3> i;
837  FTensor::Index<'J', 2> J;
838 
839  for (auto gg = 0; gg != pts.size2(); ++gg) {
840  t_gauss_pt(J) = t_ele_shape(i) * t_ele_coords(i, J);
841  t_gauss_weight = area * pts(2, gg);
842  ++t_gauss_pt;
843  ++t_gauss_weight;
844  ++t_ele_shape;
845  }
846 
847  }
848 
849  return new_pts;
850 }
851 
852 MatrixDouble
853 Tools::refineTriangleIntegrationPts(int rule, RefineTrianglesReturn refined) {
854 
855  MatrixDouble gauss_pts;
856 
857  const size_t nb_gauss_pts = QUAD_2D_TABLE[rule]->npoints;
858  gauss_pts.resize(3, nb_gauss_pts, false);
859  cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[rule]->points[1], 3,
860  &gauss_pts(0, 0), 1);
861  cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[rule]->points[2], 3,
862  &gauss_pts(1, 0), 1);
863  cblas_dcopy(nb_gauss_pts, QUAD_2D_TABLE[rule]->weights, 1, &gauss_pts(2, 0),
864  1);
865 
866  return Tools::refineTriangleIntegrationPts(gauss_pts, refined);
867 }
868 
869 } // namespace MoFEM
UBlasMatrix< double >
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::Types::VectorDouble6
VectorBoundedArray< double, 6 > VectorDouble6
Definition: Types.hpp:95
MoFEM::Tools::query_interface
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
Definition: Tools.cpp:9
MoFEM::Tools::cOre
MoFEM::Core & cOre
Definition: Tools.hpp:24
MoFEM::Tools::diffShapeFunMBEDGE
static constexpr std::array< double, 2 > diffShapeFunMBEDGE
Definition: Tools.hpp:68
MoFEM::Tools::minDistancePointFromOnSegment
static SEGMENT_MIN_DISTANCE minDistancePointFromOnSegment(const double *w_ptr, const double *v_ptr, const double *p_ptr, double *const t_ptr=nullptr)
Find closet point on the segment from the point.
Definition: Tools.cpp:446
MoFEM::Tools::SEGMENT_MIN_DISTANCE
SEGMENT_MIN_DISTANCE
Definition: Tools.hpp:525
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
EntityHandle
MoFEM::Tools::NO_SOLUTION
@ NO_SOLUTION
Definition: Tools.hpp:530
MoFEM::Tools::shapeFunMBEDGE0At00
static constexpr double shapeFunMBEDGE0At00
Definition: Tools.hpp:55
ShapeMBTET_inverse
PetscErrorCode ShapeMBTET_inverse(double *N, double *diffN, const double *elem_coords, const double *glob_coords, double *loc_coords)
calculate local coordinates for given global coordinates
Definition: fem_tools.c:335
MoFEM::Tools
Auxiliary tools.
Definition: Tools.hpp:19
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:609
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
MoFEM::Tools::checkIfPointIsInTet
static MoFEMErrorCode checkIfPointIsInTet(const double tet_coords[], const double global_coord[], const double tol, bool &result)
Check of point is in tetrahedral.
Definition: Tools.cpp:287
QUAD_2D_TABLE
static QUAD *const QUAD_2D_TABLE[]
Definition: quad.h:175
MoFEM::Tools::getTetsWithQuality
MoFEMErrorCode getTetsWithQuality(Range &out_tets, const Range &tets, Tag th=nullptr, boost::function< bool(double)> f=[](double q) -> bool { if(q<=0) return true;else return false;})
Get the Tets With Quality.
Definition: Tools.cpp:244
lapack_dgesv
static __CLPK_integer lapack_dgesv(__CLPK_integer n, __CLPK_integer nrhs, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *b, __CLPK_integer ldb)
Definition: lapack_wrap.h:176
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
FTensor::Tensor1::l2
T l2() const
Definition: Tensor1_value.hpp:84
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
nb_levels
constexpr int nb_levels
Definition: level_set.cpp:58
MoFEM::Tools::shapeFunMBTRI1At00
static constexpr double shapeFunMBTRI1At00
Definition: Tools.hpp:113
MoFEM::PetscLocalIdx_mi_tag
Definition: TagMultiIndices.hpp:45
MoFEM::Tools::diffShapeFunMBHEXAtCenter
static constexpr std::array< double, 24 > diffShapeFunMBHEXAtCenter
Definition: Tools.hpp:218
MoFEM::th
Tag th
Definition: Projection10NodeCoordsOnField.cpp:122
MoFEM::getLocalCoordinatesOnReferenceThreeNodeTriImpl
MoFEMErrorCode getLocalCoordinatesOnReferenceThreeNodeTriImpl(const T1 *elem_coords, const T2 *global_coords, const int nb_nodes, typename FTensor::promote< T1, T2 >::V *local_coords)
Definition: Tools.cpp:141
FTensor::Kronecker_Delta
Kronecker Delta class.
Definition: Kronecker_Delta.hpp:15
MoFEM::Tools::diffShapeFunMBTET
static constexpr std::array< double, 12 > diffShapeFunMBTET
Definition: Tools.hpp:271
MoFEM::Tools::shapeFunMBEDGEAt00
static constexpr std::array< double, 2 > shapeFunMBEDGEAt00
Array of shape function at zero local point on reference element.
Definition: Tools.hpp:62
J
FTensor::Index< 'J', DIM1 > J
Definition: level_set.cpp:30
MoFEM::Tools::getTriArea
double getTriArea(const EntityHandle tri) const
Get triangle area.
Definition: Tools.cpp:408
MoFEM::Tools::outerProductOfEdgeIntegrationPtsForHex
static MoFEMErrorCode outerProductOfEdgeIntegrationPtsForHex(MatrixDouble &pts, const int edge0, const int edge1, const int edge2)
Definition: Tools.cpp:659
MoFEM::Tools::RefineTrianglesReturn
std::tuple< std::vector< double >, std::vector< int >, std::vector< int > > RefineTrianglesReturn
Definition: Tools.hpp:639
FTensor::Tensor2_symmetric
Definition: Tensor2_symmetric_value.hpp: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
zeta
double zeta
Viscous hardening.
Definition: plastic.cpp:126
FTensor::Tensor2< double, 3, 3 >
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::Tools::getLocalCoordinatesOnReferenceEdgeNodeEdge
static MoFEMErrorCode getLocalCoordinatesOnReferenceEdgeNodeEdge(const double *elem_coords, const double *glob_coords, const int nb_nodes, double *local_coords)
Get the local coordinates on reference four node tet object.
Definition: Tools.cpp:203
MoFEM::Tools::minDistanceFromSegments
static SEGMENT_MIN_DISTANCE minDistanceFromSegments(const double *w_ptr, const double *v_ptr, const double *k_ptr, const double *l_ptr, double *const tvw_ptr=nullptr, double *const tlk_ptr=nullptr)
Find points on two segments in closest distance.
Definition: Tools.cpp:469
ROW
@ ROW
Definition: definitions.h:136
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
MoFEM::Tools::refineTriangleIntegrationPts
static MatrixDouble refineTriangleIntegrationPts(MatrixDouble pts, RefineTrianglesReturn refined)
generate integration points for refined triangle mesh for last level
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
QUAD_1D_TABLE
static QUAD *const QUAD_1D_TABLE[]
Definition: quad.h:164
MoFEM::Tools::getLocalCoordinatesOnReferenceThreeNodeTri
static MoFEMErrorCode getLocalCoordinatesOnReferenceThreeNodeTri(const double *elem_coords, const double *glob_coords, const int nb_nodes, double *local_coords)
Get the local coordinates on reference three node tri object.
Definition: Tools.cpp:188
FTensor::Number< 0 >
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::Tools::uniformTriangleRefineTriangles
static constexpr std::array< int, 12 > uniformTriangleRefineTriangles
Definition: Tools.hpp:629
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
a
constexpr double a
Definition: approx_sphere.cpp:30
double
RowColData
RowColData
RowColData.
Definition: definitions.h:136
MoFEM::Tools::SEGMENT_ONE_IS_POINT
@ SEGMENT_ONE_IS_POINT
Definition: Tools.hpp:527
FTensor::promote::V
T1 V
Definition: promote.hpp:17
MoFEM::Tools::findMinDistanceFromTheEdges
MoFEMErrorCode findMinDistanceFromTheEdges(const double *v_ptr, const int nb, Range edges, double *min_dist_ptr, double *o_ptr=nullptr, EntityHandle *o_segments=nullptr) const
Find minimal distance to edges.
Definition: Tools.cpp:535
COL
@ COL
Definition: definitions.h:136
eta
double eta
Definition: free_surface.cpp:170
MoFEM::Problem::getNumeredRowDofsPtr
auto & getNumeredRowDofsPtr() const
get access to numeredRowDofsPtr storing DOFs on rows
Definition: ProblemsMultiIndices.hpp:82
QUAD_::npoints
int npoints
Definition: quad.h:29
MoFEM::Tools::diffShapeFunMBQUADAtCenter
static constexpr std::array< double, 8 > diffShapeFunMBQUADAtCenter
Definition: Tools.hpp:212
MoFEM::Tools::diffShapeFunMBTRI
static constexpr std::array< double, 6 > diffShapeFunMBTRI
Definition: Tools.hpp:104
QUAD_::points
double * points
Definition: quad.h:30
MoFEM::type_from_handle
auto type_from_handle(const EntityHandle h)
get type from entity handle
Definition: Templates.hpp:1898
MoFEM::Tools::shapeFunMBTRI0At00
static constexpr double shapeFunMBTRI0At00
Definition: Tools.hpp:112
FTensor::Tensor4
Definition: Tensor4_value.hpp:18
MoFEM::Tools::getLocalCoordinatesOnReferenceFourNodeTet
static MoFEMErrorCode getLocalCoordinatesOnReferenceFourNodeTet(const double *elem_coords, const double *glob_coords, const int nb_nodes, double *local_coords)
Get the local coordinates on reference four node tet object.
Definition: Tools.cpp:88
MoFEM::Problem::getNbLocalDofsRow
DofIdx getNbLocalDofsRow() const
Definition: ProblemsMultiIndices.hpp:378
t
constexpr double t
plate stiffness
Definition: plate.cpp:58
MoFEM::Tools::SEGMENT_TWO_IS_POINT
@ SEGMENT_TWO_IS_POINT
Definition: Tools.hpp:528
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
FTensor::Index< 'i', 3 >
convert.n
n
Definition: convert.py:82
N
const int N
Definition: speed_test.cpp:3
MoFEM::UnknownInterface
base class for all interface classes
Definition: UnknownInterface.hpp:34
MoFEM::Tools::SEGMENT_TWO_AND_TWO_ARE_POINT
@ SEGMENT_TWO_AND_TWO_ARE_POINT
Definition: Tools.hpp:529
edge_coords
static const double edge_coords[6][6]
Definition: forces_and_sources_h1_continuity_check.cpp:18
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1540
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
MoFEM::Tools::minTetsQuality
MoFEMErrorCode minTetsQuality(const Range &tets, double &min_quality, Tag th=nullptr, boost::function< double(double, double)> f=[](double a, double b) -> double { return std::min(a, b);})
calculate minimal quality of tetrahedra in range
Definition: Tools.cpp:49
MoFEM::Tools::tetVolume
static double tetVolume(const double *coords)
Calculate volume of tetrahedron.
Definition: Tools.cpp:30
Range
QUAD_1D_TABLE_SIZE
#define QUAD_1D_TABLE_SIZE
Definition: quad.h:163
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
FTensor::Tensor0
Definition: Tensor0.hpp:16
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
MoFEM::Tools::refineTriangle
static RefineTrianglesReturn refineTriangle(int nb_levels)
create uniform triangle mesh of refined elements
Definition: Tools.cpp:721
MoFEM::Tools::checkVectorForNotANumber
MoFEMErrorCode checkVectorForNotANumber(const Problem *prb_ptr, const RowColData row_or_col, Vec v)
Print all DOFs for which element of vector is not a number.
Definition: Tools.cpp:307
MoFEM::Tools::writeTetsWithQuality
MoFEMErrorCode writeTetsWithQuality(const char *file_name, const char *file_type, const char *options, const Range &tets, Tag th=nullptr, boost::function< bool(double)> f=[](double q) -> bool { if(q<=0) return true;else return false;})
Write file with tetrahedral of given quality.
Definition: Tools.cpp:269
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
QUAD_::weights
double * weights
Definition: quad.h:31
MatrixBoundedArray< double, 9 >
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
ShapeDiffMBTET
PetscErrorCode ShapeDiffMBTET(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:319
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::Tools::shapeFunMBEDGE1At00
static constexpr double shapeFunMBEDGE1At00
Definition: Tools.hpp:56
MoFEM::Tools::volumeLengthQuality
static double volumeLengthQuality(const double *coords)
Calculate tetrahedron volume length quality.
Definition: Tools.cpp:15
MoFEM::Tools::getEdgeLength
static double getEdgeLength(const double *edge_coords)
Get edge length.
Definition: Tools.cpp:415
MoFEM::Problem::getNbLocalDofsCol
DofIdx getNbLocalDofsCol() const
Definition: ProblemsMultiIndices.hpp:379
MoFEM::Problem::getNumeredColDofsPtr
auto & getNumeredColDofsPtr() const
get access to numeredColDofsPtr storing DOFs on cols
Definition: ProblemsMultiIndices.hpp:87
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
MoFEM::Tools::shapeFunMBTRI2At00
static constexpr double shapeFunMBTRI2At00
Definition: Tools.hpp:114
MoFEM::Tools::outerProductOfEdgeIntegrationPtsForQuad
static MoFEMErrorCode outerProductOfEdgeIntegrationPtsForQuad(MatrixDouble &pts, const int edge0, const int edge1)
Definition: Tools.cpp:610
MoFEM::Problem
keeps basic data about problem
Definition: ProblemsMultiIndices.hpp:54
MoFEM::Tools::shapeFunMBTETAt000
static constexpr std::array< double, 4 > shapeFunMBTETAt000
Array of shape function at zero local point on reference element.
Definition: Tools.hpp:330
MoFEM::Tools::shapeFunMBTRIAt00
static constexpr std::array< double, 3 > shapeFunMBTRIAt00
Array of shape function at zero local point on reference element.
Definition: Tools.hpp:120
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MoFEM::Tools::getTriNormal
static MoFEMErrorCode getTriNormal(const double *coords, double *normal, double *d_normal=nullptr)
Get the Tri Normal objectGet triangle normal.
Definition: Tools.cpp:353
MoFEM::Tools::SOLUTION_EXIST
@ SOLUTION_EXIST
Definition: Tools.hpp:526
MoFEM::getFTensor3FromPtr< 3, 3, 3 >
FTensor::Tensor3< FTensor::PackPtr< double *, 27 >, 3, 3, 3 > getFTensor3FromPtr< 3, 3, 3 >(double *ptr)
Definition: Templates.hpp:1003
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
MOFEM_INVALID_DATA
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
tol
double tol
Definition: mesh_smoothing.cpp:26
ShapeMBTET
PetscErrorCode ShapeMBTET(double *N, const double *G_X, const double *G_Y, const double *G_Z, int DIM)
calculate shape functions
Definition: fem_tools.c:306
quad.h
Generated by Doxygen 1.8.17 and hosted at University of Glasgow