v0.15.0
Loading...
Searching...
No Matches
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
7namespace MoFEM {
8
9MoFEMErrorCode Tools::query_interface(boost::typeindex::type_index type_index,
10 UnknownInterface **iface) const {
11 *iface = const_cast<Tools *>(this);
12 return 0;
13}
14
15double 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
30double 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
48MoFEMErrorCode
49Tools::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
72constexpr double Tools::shapeFunMBEDGE0At00;
73constexpr double Tools::shapeFunMBEDGE1At00;
74constexpr std::array<double, 2> Tools::shapeFunMBEDGEAt00;
75
76constexpr std::array<double, 2> Tools::diffShapeFunMBEDGE;
77constexpr std::array<double, 6> Tools::diffShapeFunMBTRI;
78constexpr std::array<double, 12> Tools::diffShapeFunMBTET;
79constexpr double Tools::shapeFunMBTRI0At00;
80constexpr double Tools::shapeFunMBTRI1At00;
81constexpr double Tools::shapeFunMBTRI2At00;
82constexpr std::array<double, 3> Tools::shapeFunMBTRIAt00;
83
84constexpr std::array<double, 4> Tools::shapeFunMBTETAt000;
85constexpr std::array<double, 8> Tools::diffShapeFunMBQUADAtCenter;
86constexpr std::array<double, 24> Tools::diffShapeFunMBHEXAtCenter;
87
88MoFEMErrorCode 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 SETERRQ(PETSC_COMM_SELF, 1, "info == %d", info);
136
137 MoFEMFunctionReturnHot(0);
138}
139
140template <typename T1, typename T2>
141MoFEMErrorCode 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
188MoFEMErrorCode 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
195MoFEMErrorCode 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
203MoFEMErrorCode 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
244MoFEMErrorCode 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
269MoFEMErrorCode 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
287MoFEMErrorCode 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
307MoFEMErrorCode 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 SETERRQ(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
353MoFEMErrorCode 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
391MoFEMErrorCode 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
408double Tools::getTriArea(const EntityHandle tri) const {
409 if (type_from_handle(tri) != MBTRI) {
410 CHK_THROW_MESSAGE(MOFEM_INVALID_DATA, "Works only for triangle");
411 }
412 FTensor::Index<'i', 3> i;
413 FTensor::Tensor1<double, 3> t_normal;
414 CHK_THROW_MESSAGE(getTriNormal(tri, &t_normal(0)), "calculate area");
415 return sqrt(t_normal(i) * t_normal(i)) * 0.5;
416}
417
418double Tools::getEdgeLength(const double *edge_coords) {
419 FTensor::Tensor1<double, 3> t_coords_n0(edge_coords[0], edge_coords[1],
420 edge_coords[2]);
421 FTensor::Tensor1<double, 3> t_coords_n1(edge_coords[3], edge_coords[4],
422 edge_coords[5]);
423 FTensor::Index<'i', 3> i;
424 t_coords_n0(i) -= t_coords_n1(i);
425 return sqrt(t_coords_n0(i) * t_coords_n0(i));
426}
427
428double Tools::getEdgeLength(const EntityHandle edge) {
429 MoFEM::Interface &m_field = cOre;
430 moab::Interface &moab(m_field.get_moab());
431 auto get_edge_coords = [edge, &moab](double *const coords) {
432 MoFEMFunctionBegin;
433 if (type_from_handle(edge) != MBEDGE) {
434 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "Works only for edge");
435 }
436 const EntityHandle *conn;
437 int num_nodes;
438 CHKERR moab.get_connectivity(edge, conn, num_nodes, true);
439 CHKERR moab.get_coords(conn, 2, coords);
440 MoFEMFunctionReturn(0);
441 };
442 double coords[6];
443 ierr = get_edge_coords(coords);
444 CHKERRABORT(PETSC_COMM_SELF, ierr);
445 return getEdgeLength(coords);
446}
447
448Tools::SEGMENT_MIN_DISTANCE
449Tools::minDistancePointFromOnSegment(const double *w_ptr, const double *v_ptr,
450 const double *p_ptr, double *const t_ptr) {
451 FTensor::Index<'i', 3> i;
452 FTensor::Tensor1<const double *, 3> t_w(&w_ptr[0], &w_ptr[1], &w_ptr[2]);
453 FTensor::Tensor1<const double *, 3> t_v(&v_ptr[0], &v_ptr[1], &v_ptr[2]);
454 FTensor::Tensor1<const double *, 3> t_p(&p_ptr[0], &p_ptr[1], &p_ptr[2]);
455 FTensor::Tensor1<double, 3> t_vw;
456 t_vw(i) = t_v(i) - t_w(i);
457 const double dot_vw = t_vw(i) * t_vw(i);
458 if (std::abs(dot_vw) < std::numeric_limits<double>::epsilon()) {
459 if (t_ptr)
460 *t_ptr = 0;
461 return SEGMENT_ONE_IS_POINT;
462 }
463 FTensor::Tensor1<double, 3> t_pw;
464 t_pw(i) = t_p(i) - t_w(i);
465 const double t = t_pw(i) * t_vw(i) / dot_vw;
466 if (t_ptr)
467 *t_ptr = t;
468 return SOLUTION_EXIST;
469}
470
471Tools::SEGMENT_MIN_DISTANCE
472Tools::minDistanceFromSegments(const double *w_ptr, const double *v_ptr,
473 const double *k_ptr, const double *l_ptr,
474 double *const tvw_ptr, double *const tlk_ptr) {
475
476 FTensor::Index<'i', 3> i;
477 FTensor::Tensor1<const double *, 3> t_w(&w_ptr[0], &w_ptr[1], &w_ptr[2]);
478 FTensor::Tensor1<const double *, 3> t_v(&v_ptr[0], &v_ptr[1], &v_ptr[2]);
479 FTensor::Tensor1<const double *, 3> t_k(&k_ptr[0], &k_ptr[1], &k_ptr[2]);
480 FTensor::Tensor1<const double *, 3> t_l(&l_ptr[0], &l_ptr[1], &l_ptr[2]);
481
482 // First segnent is a point
483 FTensor::Tensor1<double, 3> t_vw;
484 t_vw(i) = t_v(i) - t_w(i);
485 double dot_vw = t_vw(i) * t_vw(i);
486 if (std::abs(dot_vw) < std::numeric_limits<double>::epsilon()) {
487 if (tvw_ptr)
488 *tvw_ptr = 0;
489 if (minDistancePointFromOnSegment(k_ptr, l_ptr, w_ptr, tlk_ptr) ==
490 SEGMENT_ONE_IS_POINT)
491 return SEGMENT_TWO_AND_TWO_ARE_POINT;
492 else
493 return SEGMENT_ONE_IS_POINT;
494 }
495
496 // Second segment is a point
497 FTensor::Tensor1<double, 3> t_lk;
498 t_lk(i) = t_l(i) - t_k(i);
499 double dot_lk = t_lk(i) * t_lk(i);
500 if (std::abs(dot_lk) < std::numeric_limits<double>::epsilon()) {
501 if (tlk_ptr)
502 *tlk_ptr = 0;
503 if (minDistancePointFromOnSegment(w_ptr, v_ptr, k_ptr, tvw_ptr) ==
504 SEGMENT_ONE_IS_POINT)
505 return SEGMENT_TWO_AND_TWO_ARE_POINT;
506 else
507 return SEGMENT_TWO_IS_POINT;
508 }
509
510 const double a = t_vw(i) * t_vw(i);
511 const double b = -t_vw(i) * t_lk(i);
512 const double c = t_lk(i) * t_lk(i);
513
514 const double det = a * c - b * b;
515 if (std::abs(det) < std::numeric_limits<double>::epsilon()) {
516
517 return NO_SOLUTION;
518
519 } else {
520
521 FTensor::Tensor1<double, 3> t_wk;
522 t_wk(i) = t_w(i) - t_k(i);
523
524 const double ft0 = t_vw(i) * t_wk(i);
525 const double ft1 = -t_lk(i) * t_wk(i);
526 const double t0 = (ft1 * b - ft0 * c) / det;
527 const double t1 = (ft0 * b - ft1 * a) / det;
528
529 if (tvw_ptr)
530 *tvw_ptr = t0;
531 if (tlk_ptr)
532 *tlk_ptr = t1;
533
534 return SOLUTION_EXIST;
535 }
536}
537
538MoFEMErrorCode Tools::findMinDistanceFromTheEdges(
539 const double *v_ptr, const int nb, Range edges, double *min_dist_ptr,
540 double *o_ptr, EntityHandle *o_segments) const {
541 MoFEM::Interface &m_field = cOre;
542 moab::Interface &moab(m_field.get_moab());
543 MoFEMFunctionBegin;
544
545 FTensor::Index<'i', 3> i;
546
547 auto get_point = [i](auto &t_w, auto &t_delta, auto t) {
548 FTensor::Tensor1<double, 3> t_p;
549 t = std::max(0., std::min(1., t));
550 t_p(i) = t_w(i) + t * t_delta(i);
551 return t_p;
552 };
553
554 auto get_distance = [i](auto &t_p, auto &t_n) {
555 FTensor::Tensor1<double, 3> t_dist_vector;
556 t_dist_vector(i) = t_p(i) - t_n(i);
557 return sqrt(t_dist_vector(i) * t_dist_vector(i));
558 };
559
560 for (auto e : edges) {
561
562 int num_nodes;
563 const EntityHandle *conn_fixed;
564 CHKERR moab.get_connectivity(e, conn_fixed, num_nodes, true);
565 VectorDouble6 coords_fixed(6);
566 CHKERR moab.get_coords(conn_fixed, num_nodes, &coords_fixed[0]);
567 FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_f0(
568 &coords_fixed[0], &coords_fixed[1], &coords_fixed[2]);
569 FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_f1(
570 &coords_fixed[3], &coords_fixed[4], &coords_fixed[5]);
571
572 FTensor::Tensor1<double, 3> t_edge_delta;
573 t_edge_delta(i) = t_f1(i) - t_f0(i);
574
575 FTensor::Tensor1<FTensor::PackPtr<const double *, 3>, 3> t_n(
576 v_ptr, v_ptr + 1, v_ptr + 2);
577 FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_min_coords(
578 o_ptr, o_ptr + 1, o_ptr + 2);
579 FTensor::Tensor0<FTensor::PackPtr<double *, 1>> t_min_dist(min_dist_ptr);
580
581 EntityHandle *colsest_segment_it = nullptr;
582 if (o_segments)
583 colsest_segment_it = o_segments;
584
585 for (int n = 0; n != nb; ++n) {
586
587 double t;
588 if (Tools::minDistancePointFromOnSegment(&t_f0(0), &t_f1(0), &t_n(0),
589 &t) == Tools::SOLUTION_EXIST) {
590 auto t_p = get_point(t_f0, t_edge_delta, t);
591 auto dist_n = get_distance(t_p, t_n);
592 if (dist_n < t_min_dist || t_min_dist < 0) {
593 t_min_dist = dist_n;
594 if (o_ptr)
595 t_min_coords(i) = t_p(i);
596 if (o_segments)
597 *colsest_segment_it = e;
598 }
599 }
600
601 ++t_n;
602 ++t_min_dist;
603 if (o_ptr)
604 ++t_min_coords;
605 if (o_segments)
606 ++colsest_segment_it;
607 }
608 }
609
610 MoFEMFunctionReturn(0);
611}
612
613MoFEMErrorCode Tools::outerProductOfEdgeIntegrationPtsForQuad(
614 MatrixDouble &gauss_pts, const int rule_ksi, const int rule_eta) {
615 MoFEMFunctionBegin;
616
617 auto check_rule_edge = [](int rule) {
618 MoFEMFunctionBeginHot;
619 if (rule < 0) {
620 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
621 "Wrong integration rule: %d", rule);
622 }
623 if (rule > QUAD_1D_TABLE_SIZE) {
624 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
625 "rule > quadrature order %d < %d", rule, QUAD_1D_TABLE_SIZE);
626 }
627 if (QUAD_1D_TABLE[rule]->dim != 1) {
628 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong dimension");
629 }
630 if (QUAD_1D_TABLE[rule]->order < rule) {
631 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
632 "wrong order %d != %d", QUAD_1D_TABLE[rule]->order, rule);
633 }
634 MoFEMFunctionReturnHot(0);
635 };
636
637 CHKERR check_rule_edge(rule_ksi);
638 CHKERR check_rule_edge(rule_eta);
639
640 const int nb_gauss_pts_ksi = QUAD_1D_TABLE[rule_ksi]->npoints;
641 const int nb_gauss_pts_eta = QUAD_1D_TABLE[rule_eta]->npoints;
642 gauss_pts.resize(3, nb_gauss_pts_ksi * nb_gauss_pts_eta, false);
643
644 int gg = 0;
645 for (size_t i = 0; i != nb_gauss_pts_ksi; ++i) {
646 const double wi = QUAD_1D_TABLE[rule_ksi]->weights[i];
647 const double ksi = (QUAD_1D_TABLE[rule_ksi]->points[2 * i + 1]);
648 for (size_t j = 0; j != nb_gauss_pts_eta; ++j, ++gg) {
649 const double wk = wi * QUAD_1D_TABLE[rule_eta]->weights[j];
650 const double eta = QUAD_1D_TABLE[rule_eta]->points[2 * j + 1];
651 gauss_pts(0, gg) = ksi;
652 gauss_pts(1, gg) = eta;
653 gauss_pts(2, gg) = wk;
654 }
655 }
656 if (gg != gauss_pts.size2())
657 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong size of matrix");
658
659 MoFEMFunctionReturn(0);
660}
661
662MoFEMErrorCode Tools::outerProductOfEdgeIntegrationPtsForHex(
663 MatrixDouble &gauss_pts, const int rule_ksi, const int rule_eta,
664 const int rule_zeta) {
665 MoFEMFunctionBegin;
666
667 auto check_rule_edge = [](int rule) {
668 MoFEMFunctionBeginHot;
669 if (rule < 0) {
670 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
671 "Wrong integration rule: %d", rule);
672 }
673 if (rule > QUAD_1D_TABLE_SIZE) {
674 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
675 "rule > quadrature order %d < %d", rule, QUAD_1D_TABLE_SIZE);
676 }
677 if (QUAD_1D_TABLE[rule]->dim != 1) {
678 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong dimension");
679 }
680 if (QUAD_1D_TABLE[rule]->order < rule) {
681 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
682 "wrong order %d != %d", QUAD_1D_TABLE[rule]->order, rule);
683 }
684 MoFEMFunctionReturnHot(0);
685 };
686
687 CHKERR check_rule_edge(rule_ksi);
688 CHKERR check_rule_edge(rule_eta);
689 CHKERR check_rule_edge(rule_zeta);
690
691 const int nb_gauss_pts_ksi = QUAD_1D_TABLE[rule_ksi]->npoints;
692 const int nb_gauss_pts_eta = QUAD_1D_TABLE[rule_eta]->npoints;
693 const int nb_gauss_pts_zeta = QUAD_1D_TABLE[rule_zeta]->npoints;
694 gauss_pts.resize(4, nb_gauss_pts_ksi * nb_gauss_pts_eta * nb_gauss_pts_zeta,
695 false);
696
697 int gg = 0;
698 for (size_t i = 0; i != nb_gauss_pts_ksi; ++i) {
699 const double wi = QUAD_1D_TABLE[rule_ksi]->weights[i];
700 const double ksi = QUAD_1D_TABLE[rule_ksi]->points[2 * i + 1];
701 for (size_t j = 0; j != nb_gauss_pts_eta; ++j) {
702 const double wj = wi * QUAD_1D_TABLE[rule_eta]->weights[j];
703 const double eta = QUAD_1D_TABLE[rule_eta]->points[2 * j + 1];
704 for (size_t k = 0; k != nb_gauss_pts_zeta; ++k, ++gg) {
705 const double wk = wj * QUAD_1D_TABLE[rule_zeta]->weights[k];
706 const double zeta = QUAD_1D_TABLE[rule_zeta]->points[2 * k + 1];
707 gauss_pts(0, gg) = ksi;
708 gauss_pts(1, gg) = eta;
709 gauss_pts(2, gg) = zeta;
710 gauss_pts(3, gg) = wk;
711 }
712 }
713 }
714
715 if (gg != gauss_pts.size2())
716 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong size of matrix");
717
718 MoFEMFunctionReturn(0);
719}
720
721constexpr std::array<int, 12> Tools::uniformTriangleRefineTriangles;
722
723std::tuple<std::vector<double>, std::vector<int>, std::vector<int>>
724Tools::refineTriangle(int nb_levels) {
725
726 std::vector<int> triangles{0, 1, 2, 3, 4, 5};
727 std::vector<double> nodes{
728
729 0., 0., // 0
730 1., 0., // 1
731 0., 1., // 2
732 0.5, 0., // 3
733 0.5, 0.5, // 4
734 0., 0.5 // 5
735
736 };
737 std::map<std::pair<int, int>, int> edges{
738 {{0, 1}, 3}, {{1, 2}, 4}, {{0, 2}, 5}};
739
740 auto add_edge = [&](auto a, auto b) {
741 if (a > b) {
742 std::swap(a, b);
743 }
744 auto it = edges.find(std::make_pair(a, b));
745 if (it == edges.end()) {
746 int e = 3 + edges.size();
747 edges[std::make_pair(a, b)] = e;
748 for (auto n : {0, 1}) {
749 nodes.push_back((nodes[2 * a + n] + nodes[2 * b + n]) / 2);
750 }
751#ifndef NDEBUG
752 if (e != nodes.size() / 2 - 1)
753 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "wrong node/edge index");
754#endif
755 return e;
756 }
757 return it->second;
758 };
759
760 auto add_triangle = [&](auto t) {
761 for (auto tt : {0, 1, 2, 3}) {
762 auto last = triangles.size() / 6;
763 for (auto n : {0, 1, 2}) {
764 // add triangle nodes
765 triangles.push_back(
766 triangles[6 * t + uniformTriangleRefineTriangles[3 * tt + n]]);
767 }
768 // add triangle edges
769 auto cycle = std::array<int, 4>{0, 1, 2, 0};
770 for (auto e : {0, 1, 2}) {
771 triangles.push_back(add_edge(triangles[6 * last + cycle[e]],
772 triangles[6 * last + cycle[e + 1]]));
773 }
774 }
775 };
776
777 std::vector<int> level_index{0, 1};
778 auto l = 0;
779 for (; l != nb_levels; ++l) {
780 auto first_tet = level_index[l];
781 auto nb_last_level_test = level_index.back() - level_index[l];
782 for (auto t = first_tet; t != (first_tet + nb_last_level_test); ++t) {
783 add_triangle(t);
784 }
785 level_index.push_back(triangles.size() / 6);
786 }
787
788 return std::make_tuple(nodes, triangles, level_index);
789}
790
791MatrixDouble Tools::refineTriangleIntegrationPts(
792 MatrixDouble pts,
793 std::tuple<std::vector<double>, std::vector<int>, std::vector<int>>
794 refined) {
795
796 auto [nodes, triangles, level_index] = refined;
797
798 auto get_coords = [&](auto t) {
799 auto [nodes, triangles, level_index] = refined;
800 std::array<double, 9> ele_coords;
801 for (auto n : {0, 1, 2}) {
802 for (auto i : {0, 1}) {
803 ele_coords[3 * n + i] = nodes[2 * triangles[6 * t + n] + i];
804 }
805 ele_coords[3 * n + 2] = 0;
806 }
807 return ele_coords;
808 };
809
810 auto get_normal = [](auto &ele_coords) {
811 FTensor::Tensor1<double, 3> t_normal;
812 Tools::getTriNormal(ele_coords.data(), &t_normal(0));
813 return t_normal;
814 };
815
816 std::vector<double> ele_shape(3 * pts.size2());
817 shapeFunMBTRI<1>(&*ele_shape.begin(), &pts(0, 0), &pts(1, 0), pts.size2());
818
819 int nb_elems = level_index.back() - level_index[level_index.size() - 2];
820 MatrixDouble new_pts(3, pts.size2() * nb_elems);
821 FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 2> t_gauss_pt{&new_pts(0, 0),
822 &new_pts(1, 0)};
823 FTensor::Tensor0<FTensor::PackPtr<double *, 1>> t_gauss_weight{
824 &new_pts(2, 0)};
825
826 for (auto t = level_index[level_index.size() - 2]; t != level_index.back();
827 ++t) {
828
829 auto ele_coords = get_coords(t);
830 auto t_normal = get_normal(ele_coords);
831 auto area = t_normal.l2();
832
833 FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_ele_shape{
834 &ele_shape[0], &ele_shape[1], &ele_shape[2]};
835 FTensor::Tensor2<double, 3, 2> t_ele_coords{ele_coords[0], ele_coords[1],
836 ele_coords[3], ele_coords[4],
837 ele_coords[6], ele_coords[7]};
838
839 FTensor::Index<'i', 3> i;
840 FTensor::Index<'J', 2> J;
841
842 for (auto gg = 0; gg != pts.size2(); ++gg) {
843 t_gauss_pt(J) = t_ele_shape(i) * t_ele_coords(i, J);
844 t_gauss_weight = area * pts(2, gg);
845 ++t_gauss_pt;
846 ++t_gauss_weight;
847 ++t_ele_shape;
848 }
849
850 }
851
852 return new_pts;
853}
854
855MatrixDouble
856Tools::refineTriangleIntegrationPts(int rule, RefineTrianglesReturn refined) {
857
858 MatrixDouble gauss_pts;
859
860 const size_t nb_gauss_pts = QUAD_2D_TABLE[rule]->npoints;
861 gauss_pts.resize(3, nb_gauss_pts, false);
862 cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[rule]->points[1], 3,
863 &gauss_pts(0, 0), 1);
864 cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[rule]->points[2], 3,
865 &gauss_pts(1, 0), 1);
866 cblas_dcopy(nb_gauss_pts, QUAD_2D_TABLE[rule]->weights, 1, &gauss_pts(2, 0),
867 1);
868
869 return Tools::refineTriangleIntegrationPts(gauss_pts, refined);
870}
871
872} // namespace MoFEM
a
constexpr double a
Definition approx_sphere.cpp:30
FTensor::Kronecker_Delta
Kronecker Delta class.
Definition Kronecker_Delta.hpp:15
FTensor::Tensor2_symmetric
Definition Tensor2_symmetric_value.hpp:14
RowColData
RowColData
RowColData.
Definition definitions.h:136
COL
@ COL
Definition definitions.h:136
ROW
@ ROW
Definition definitions.h:136
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition definitions.h:612
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:463
L2
@ L2
field with C-1 continuity
Definition definitions.h:88
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
MOFEM_INVALID_DATA
@ MOFEM_INVALID_DATA
Definition definitions.h:36
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:432
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:456
order
constexpr int order
Definition dg_projection.cpp:18
ShapeDiffMBTET
PetscErrorCode ShapeDiffMBTET(double *diffN)
calculate derivatives of shape functions
Definition fem_tools.c:319
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
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
edge_coords
static const double edge_coords[6][6]
Definition forces_and_sources_h1_continuity_check.cpp:18
t_kd
constexpr auto t_kd
Definition free_surface.cpp:139
eta
double eta
Definition free_surface.cpp:172
i
FTensor::Index< 'i', SPACE_DIM > i
Definition hcurl_divergence_operator_2d.cpp:27
c
const double c
speed of light (cm/ns)
Definition initial_diffusion.cpp:39
v
const double v
phase velocity of light in medium (cm/ns)
Definition initial_diffusion.cpp:40
n
const double n
refractive index of diffusive medium
Definition initial_diffusion.cpp:38
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:188
nb_levels
constexpr int nb_levels
Definition level_set.cpp:58
J
FTensor::Index< 'J', DIM1 > J
Definition level_set.cpp:30
l
FTensor::Index< 'l', 3 > l
Definition matrix_function.cpp:21
j
FTensor::Index< 'j', 3 > j
Definition matrix_function.cpp:19
k
FTensor::Index< 'k', 3 > k
Definition matrix_function.cpp:20
tol
double tol
Definition mesh_smoothing.cpp:26
FTensor::levi_civita
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
Definition Levi_Civita.hpp:617
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition Exceptions.hpp:76
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition Exceptions.hpp:56
MoFEM::Types::VectorDouble6
VectorBoundedArray< double, 6 > VectorDouble6
Definition Types.hpp:95
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition Types.hpp:77
MoFEM
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
MoFEM::type_from_handle
auto type_from_handle(const EntityHandle h)
get type from entity handle
Definition Templates.hpp:1932
MoFEM::getFTensor2FromPtr
auto getFTensor2FromPtr(double *ptr)
Make Tensor2 from pointer.
Definition Templates.hpp:925
MoFEM::getFTensor1FromPtr
FTensor::Tensor1< FTensor::PackPtr< double *, S >, DIM > getFTensor1FromPtr(double *ptr)
Make Tensor1 from pointer.
Definition Templates.hpp:838
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
MoFEM::getFTensor3FromPtr< 3, 3, 3 >
FTensor::Tensor3< FTensor::PackPtr< double *, 27 >, 3, 3, 3 > getFTensor3FromPtr< 3, 3, 3 >(double *ptr)
Definition Templates.hpp:1003
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition Templates.hpp:1574
zeta
double zeta
Viscous hardening.
Definition plastic.cpp:130
t
constexpr double t
plate stiffness
Definition plate.cpp:58
QUAD_2D_TABLE
static QUAD *const QUAD_2D_TABLE[]
Definition quad.h:175
QUAD_1D_TABLE
static QUAD *const QUAD_1D_TABLE[]
Definition quad.h:164
QUAD_1D_TABLE_SIZE
#define QUAD_1D_TABLE_SIZE
Definition quad.h:163
m
FTensor::Index< 'm', 3 > m
Definition shallow_wave.cpp:80
N
const int N
Definition speed_test.cpp:3
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition DeprecatedCoreInterface.hpp:16
MoFEM::PetscLocalIdx_mi_tag
Definition TagMultiIndices.hpp:45
MoFEM::Problem
keeps basic data about problem
Definition ProblemsMultiIndices.hpp:54
MoFEM::Problem::getNbLocalDofsRow
DofIdx getNbLocalDofsRow() const
Definition ProblemsMultiIndices.hpp:378
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
MoFEM::Problem::getNumeredRowDofsPtr
auto & getNumeredRowDofsPtr() const
get access to numeredRowDofsPtr storing DOFs on rows
Definition ProblemsMultiIndices.hpp:82
MoFEM::Tools
Auxiliary tools.
Definition Tools.hpp:19
MoFEM::Tools::shapeFunMBEDGE1At00
static constexpr double shapeFunMBEDGE1At00
Definition Tools.hpp:56
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:472
MoFEM::Tools::refineTriangle
static RefineTrianglesReturn refineTriangle(int nb_levels)
create uniform triangle mesh of refined elements
Definition Tools.cpp:724
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::Tools::outerProductOfEdgeIntegrationPtsForQuad
static MoFEMErrorCode outerProductOfEdgeIntegrationPtsForQuad(MatrixDouble &pts, const int edge0, const int edge1)
Definition Tools.cpp:613
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:449
MoFEM::Tools::shapeFunMBTRI1At00
static constexpr double shapeFunMBTRI1At00
Definition Tools.hpp:113
MoFEM::Tools::RefineTrianglesReturn
std::tuple< std::vector< double >, std::vector< int >, std::vector< int > > RefineTrianglesReturn
Definition Tools.hpp:647
MoFEM::Tools::shapeFunMBTRI0At00
static constexpr double shapeFunMBTRI0At00
Definition Tools.hpp:112
MoFEM::Tools::diffShapeFunMBQUADAtCenter
static constexpr std::array< double, 8 > diffShapeFunMBQUADAtCenter
Definition Tools.hpp:212
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:538
MoFEM::Tools::diffShapeFunMBTRI
static constexpr std::array< double, 6 > diffShapeFunMBTRI
Definition Tools.hpp:104
MoFEM::Tools::SEGMENT_TWO_AND_TWO_ARE_POINT
@ SEGMENT_TWO_AND_TWO_ARE_POINT
Definition Tools.hpp:538
MoFEM::Tools::SEGMENT_ONE_IS_POINT
@ SEGMENT_ONE_IS_POINT
Definition Tools.hpp:536
MoFEM::Tools::SEGMENT_TWO_IS_POINT
@ SEGMENT_TWO_IS_POINT
Definition Tools.hpp:537
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::shapeFunMBEDGEAt00
static constexpr std::array< double, 2 > shapeFunMBEDGEAt00
Array of shape function at zero local point on reference element.
Definition Tools.hpp:62
MoFEM::Tools::shapeFunMBEDGE0At00
static constexpr double shapeFunMBEDGE0At00
Definition Tools.hpp:55
MoFEM::Tools::refineTriangleIntegrationPts
static MatrixDouble refineTriangleIntegrationPts(MatrixDouble pts, RefineTrianglesReturn refined)
generate integration points for refined triangle mesh for last level
Definition Tools.cpp:791
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::uniformTriangleRefineTriangles
static constexpr std::array< int, 12 > uniformTriangleRefineTriangles
Definition Tools.hpp:638
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
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
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
MoFEM::Tools::outerProductOfEdgeIntegrationPtsForHex
static MoFEMErrorCode outerProductOfEdgeIntegrationPtsForHex(MatrixDouble &pts, const int edge0, const int edge1, const int edge2)
Definition Tools.cpp:662
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::shapeFunMBTRI2At00
static constexpr double shapeFunMBTRI2At00
Definition Tools.hpp:114
MoFEM::Tools::diffShapeFunMBHEXAtCenter
static constexpr std::array< double, 24 > diffShapeFunMBHEXAtCenter
Definition Tools.hpp:218
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::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
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::diffShapeFunMBEDGE
static constexpr std::array< double, 2 > diffShapeFunMBEDGE
Definition Tools.hpp:68
MoFEM::Tools::tetVolume
static double tetVolume(const double *coords)
Calculate volume of tetrahedron.
Definition Tools.cpp:30
MoFEM::Tools::volumeLengthQuality
static double volumeLengthQuality(const double *coords)
Calculate tetrahedron volume length quality.
Definition Tools.cpp:15
MoFEM::Tools::getTriArea
double getTriArea(const EntityHandle tri) const
Get triangle area.
Definition Tools.cpp:408
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
MoFEM::Tools::query_interface
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
Definition Tools.cpp:9
MoFEM::Tools::diffShapeFunMBTET
static constexpr std::array< double, 12 > diffShapeFunMBTET
Definition Tools.hpp:271
MoFEM::Tools::getEdgeLength
static double getEdgeLength(const double *edge_coords)
Get edge length.
Definition Tools.cpp:418
MoFEM::Tools::cOre
MoFEM::Core & cOre
Definition Tools.hpp:24
MoFEM::Tools::shapeFunMBTRI
static MoFEMErrorCode shapeFunMBTRI(double *shape, const double *ksi, const double *eta, const int nb)
Calculate shape functions on triangle.
Definition Tools.hpp:716
MoFEM::UnknownInterface
base class for all interface classes
Definition UnknownInterface.hpp:34
QUAD_::points
double * points
Definition quad.h:30
QUAD_::weights
double * weights
Definition quad.h:31
QUAD_::npoints
int npoints
Definition quad.h:29
Generated by Doxygen 1.12.0 and hosted at University of Glasgow