v0.13.2
Loading...
Searching...
No Matches
Tools.cpp
Go to the documentation of this file.
1/** \file Tools.cpp
2 * \brief Auxiliary tools
3 */
4
5
6
7#include <phg-quadrule/quad.h>
8
9namespace MoFEM {
10
11MoFEMErrorCode Tools::query_interface(boost::typeindex::type_index type_index,
12 UnknownInterface **iface) const {
13 *iface = const_cast<Tools *>(this);
14 return 0;
15}
16
17double Tools::volumeLengthQuality(const double *coords) {
18 double lrms = 0;
19 for (int dd = 0; dd != 3; dd++) {
20 lrms += pow(coords[0 * 3 + dd] - coords[1 * 3 + dd], 2) +
21 pow(coords[0 * 3 + dd] - coords[2 * 3 + dd], 2) +
22 pow(coords[0 * 3 + dd] - coords[3 * 3 + dd], 2) +
23 pow(coords[1 * 3 + dd] - coords[2 * 3 + dd], 2) +
24 pow(coords[1 * 3 + dd] - coords[3 * 3 + dd], 2) +
25 pow(coords[2 * 3 + dd] - coords[3 * 3 + dd], 2);
26 }
27 lrms = sqrt((1. / 6.) * lrms);
28 double volume = tetVolume(coords);
29 return 6. * sqrt(2.) * volume / pow(lrms, 3);
30}
31
32double Tools::tetVolume(const double *coords) {
33 double diff_n[12];
34 ShapeDiffMBTET(diff_n);
35 FTensor::Tensor1<double *, 3> t_diff_n(&diff_n[0], &diff_n[1], &diff_n[2], 3);
36 FTensor::Tensor1<const double *, 3> t_coords(&coords[0], &coords[1],
37 &coords[2], 3);
38 FTensor::Tensor2<double, 3, 3> jac;
39 FTensor::Index<'i', 3> i;
40 FTensor::Index<'j', 3> j;
41 jac(i, j) = 0;
42 for (int nn = 0; nn != 4; nn++) {
43 jac(i, j) += t_coords(i) * t_diff_n(j);
44 ++t_coords;
45 ++t_diff_n;
46 }
47 return determinantTensor3by3(jac) / 6.;
48}
49
50MoFEMErrorCode
51Tools::minTetsQuality(const Range &tets, double &min_quality, Tag th,
52 boost::function<double(double, double)> f) {
53 MoFEM::Interface &m_field = cOre;
54 moab::Interface &moab(m_field.get_moab());
55 MoFEMFunctionBegin;
56 const EntityHandle *conn;
57 int num_nodes;
58 double coords[12];
59 for (auto tet : tets) {
60 CHKERR m_field.get_moab().get_connectivity(tet, conn, num_nodes, true);
61 if (th) {
62 CHKERR moab.tag_get_data(th, conn, num_nodes, coords);
63 } else {
64 CHKERR moab.get_coords(conn, num_nodes, coords);
65 }
66 double q = Tools::volumeLengthQuality(coords);
67 if (!std::isnormal(q))
68 q = -2;
69 min_quality = f(q, min_quality);
70 }
71 MoFEMFunctionReturn(0);
72}
73
74constexpr double Tools::shapeFunMBEDGE0At00;
75constexpr double Tools::shapeFunMBEDGE1At00;
76constexpr std::array<double, 2> Tools::shapeFunMBEDGEAt00;
77
78constexpr std::array<double, 2> Tools::diffShapeFunMBEDGE;
79constexpr std::array<double, 6> Tools::diffShapeFunMBTRI;
80constexpr std::array<double, 12> Tools::diffShapeFunMBTET;
81constexpr double Tools::shapeFunMBTRI0At00;
82constexpr double Tools::shapeFunMBTRI1At00;
83constexpr double Tools::shapeFunMBTRI2At00;
84constexpr std::array<double, 3> Tools::shapeFunMBTRIAt00;
85
86constexpr std::array<double, 4> Tools::shapeFunMBTETAt000;
87constexpr std::array<double, 8> Tools::diffShapeFunMBQUADAtCenter;
88constexpr std::array<double, 24> Tools::diffShapeFunMBHEXAtCenter;
89
90MoFEMErrorCode Tools::getLocalCoordinatesOnReferenceFourNodeTet(
91 const double *elem_coords, const double *global_coords, const int nb_nodes,
92 double *local_coords) {
93 FTensor::Index<'i', 4> i;
94 MoFEMFunctionBeginHot;
95
96 FTensor::Tensor1<FTensor::PackPtr<const double *, 1>, 4> t_elem_coords = {
97 &elem_coords[0], &elem_coords[3], &elem_coords[6], &elem_coords[9]};
98
99 FTensor::Tensor1<const double, 4> t_n = {
100 shapeFunMBTETAt000[0], shapeFunMBTETAt000[1], shapeFunMBTETAt000[2],
101 shapeFunMBTETAt000[3]};
102 FTensor::Tensor1<double, 3> t_coords_at_0;
103
104 // Build matrix and get coordinates of zero point
105 // ii - global coordinates
106 // jj - local derivatives
107 MatrixDouble3by3 a(3, 3);
108 for (auto ii : {0, 1, 2}) {
109 FTensor::Tensor1<FTensor::PackPtr<const double *, 1>, 4> t_diff(
110 &diffShapeFunMBTET[0], &diffShapeFunMBTET[3], &diffShapeFunMBTET[6],
111 &diffShapeFunMBTET[9]);
112 for (auto jj : {0, 1, 2}) {
113 a(jj, ii) = t_diff(i) * t_elem_coords(i);
114 ++t_diff;
115 }
116 t_coords_at_0(ii) = t_n(i) * t_elem_coords(i);
117 ++t_elem_coords;
118 }
119
120 FTensor::Tensor1<FTensor::PackPtr<const double *, 3>, 3> t_global_coords = {
121 &global_coords[0], &global_coords[1], &global_coords[2]};
122 FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_local_coords = {
123 &local_coords[0], &local_coords[1], &local_coords[2]};
124
125 // Calculate right hand side
126 FTensor::Index<'j', 3> j;
127 for (int ii = 0; ii != nb_nodes; ++ii) {
128 t_local_coords(j) = t_global_coords(j) - t_coords_at_0(j);
129 ++t_local_coords;
130 ++t_global_coords;
131 }
132
133 // Solve problem
134 int IPIV[3];
135 int info = lapack_dgesv(3, nb_nodes, &a(0, 0), 3, IPIV, local_coords, 3);
136 if (info != 0)
137 SETERRQ1(PETSC_COMM_SELF, 1, "info == %d", info);
138
139 MoFEMFunctionReturnHot(0);
140}
141
142MoFEMErrorCode Tools::getLocalCoordinatesOnReferenceTriNodeTri(
143 const double *elem_coords, const double *global_coords, const int nb_nodes,
144 double *local_coords) {
145
146 FTensor::Index<'i', 3> i3;
147 FTensor::Index<'j', 3> j3;
148 FTensor::Index<'K', 2> K2;
149 FTensor::Index<'L', 2> L2;
150 MoFEMFunctionBeginHot;
151
152 FTensor::Tensor1<FTensor::PackPtr<const double *, 1>, 3> t_elem_coords = {
153 &elem_coords[0], &elem_coords[3], &elem_coords[6]};
154
155 FTensor::Tensor1<const double, 3> t_n = {
156 shapeFunMBTRIAt00[0], shapeFunMBTRIAt00[1], shapeFunMBTRIAt00[2]};
157 FTensor::Tensor1<double, 3> t_coords_at_0;
158
159 // Build matrix and get coordinates of zero point
160 // ii - global coordinates
161 // jj - local derivatives
162 FTensor::Tensor2<double, 2, 3> t_a;
163 for (auto ii : {0, 1, 2}) {
164 FTensor::Tensor1<FTensor::PackPtr<const double *, 1>, 3> t_diff(
165 &diffShapeFunMBTRI[0], &diffShapeFunMBTRI[2], &diffShapeFunMBTRI[4]);
166 for (auto jj : {0, 1}) {
167 t_a(jj, ii) = t_diff(i3) * t_elem_coords(i3);
168 ++t_diff;
169 }
170 t_coords_at_0(ii) = t_n(i3) * t_elem_coords(i3);
171 ++t_elem_coords;
172 }
173
174 FTensor::Tensor1<FTensor::PackPtr<const double *, 3>, 3> t_global_coords = {
175 &global_coords[0], &global_coords[1], &global_coords[2]};
176 FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2> t_local_coords = {
177 &local_coords[0], &local_coords[1]};
178
179 FTensor::Tensor2_symmetric<double, 2> t_b, t_inv_b;
180 t_b(K2, L2) = t_a(K2, j3) ^ t_a(L2, j3);
181
182 // Solve problem
183 const auto inv_det = 1. / (t_b(0, 0) * t_b(1, 1) - t_b(0, 1) * t_b(1, 0));
184 t_inv_b(0, 0) = t_b(1, 1) * inv_det;
185 t_inv_b(0, 1) = -t_b(0, 1) * inv_det;
186 t_inv_b(1, 1) = t_b(0, 0) * inv_det;
187
188 // Calculate right hand side
189 for (int ii = 0; ii != nb_nodes; ++ii) {
190 t_local_coords(L2) =
191 t_inv_b(L2, K2) *
192 (t_a(K2, j3) * (t_global_coords(j3) - t_coords_at_0(j3)));
193 ++t_local_coords;
194 ++t_global_coords;
195 }
196
197 MoFEMFunctionReturnHot(0);
198}
199
200MoFEMErrorCode Tools::getLocalCoordinatesOnReferenceEdgeNodeEdge(
201 const double *elem_coords, const double *global_coords, const int nb_nodes,
202 double *local_coords) {
203
204 FTensor::Index<'i', 3> i3;
205 MoFEMFunctionBeginHot;
206
207 FTensor::Tensor1<FTensor::PackPtr<const double *, 1>, 3> t_elem_coords = {
208 &elem_coords[0], &elem_coords[3], &elem_coords[6]};
209
210 FTensor::Tensor1<const double, 2> t_n = {shapeFunMBEDGEAt00[0],
211 shapeFunMBEDGEAt00[1]};
212 FTensor::Tensor1<double, 3> t_coords_at_0;
213 // Build matrix and get coordinates of zero point
214 // ii - global coordinates
215 FTensor::Tensor1<double, 3> t_a;
216 for (auto ii : {0, 1, 2}) {
217 t_a(ii) = diffShapeFunMBEDGE[0] * t_elem_coords(0) +
218 diffShapeFunMBEDGE[1] * t_elem_coords(1);
219 t_coords_at_0(ii) = shapeFunMBEDGEAt00[0] * t_elem_coords(0) +
220 shapeFunMBEDGEAt00[1] * t_elem_coords(1);
221 ++t_elem_coords;
222 }
223
224 FTensor::Tensor1<FTensor::PackPtr<const double *, 3>, 3> t_global_coords = {
225 &global_coords[0], &global_coords[1], &global_coords[2]};
226 FTensor::Tensor0<FTensor::PackPtr<double *, 1>> t_local_coords = {
227 &local_coords[0]};
228
229 const double b = t_a(i3) * t_a(i3);
230 const double inv_b = 1 / b;
231
232 // Calculate right hand side
233 for (int ii = 0; ii != nb_nodes; ++ii) {
234 t_local_coords =
235 inv_b * (t_a(i3) * (t_global_coords(i3) - t_coords_at_0(i3)));
236 ++t_local_coords;
237 ++t_global_coords;
238 }
239
240 MoFEMFunctionReturnHot(0);
241}
242
243MoFEMErrorCode Tools::getTetsWithQuality(Range &out_tets, const Range &tets,
244 Tag th,
245 boost::function<bool(double)> f) {
246 MoFEM::Interface &m_field = cOre;
247 moab::Interface &moab(m_field.get_moab());
248 MoFEMFunctionBegin;
249 Range to_write;
250 const EntityHandle *conn;
251 int num_nodes;
252 double coords[12];
253 for (auto tet : tets) {
254 CHKERR m_field.get_moab().get_connectivity(tet, conn, num_nodes, true);
255 if (th) {
256 CHKERR moab.tag_get_data(th, conn, num_nodes, coords);
257 } else {
258 CHKERR moab.get_coords(conn, num_nodes, coords);
259 }
260 double q = Tools::volumeLengthQuality(coords);
261 if (f(q)) {
262 out_tets.insert(tet);
263 }
264 }
265 MoFEMFunctionReturn(0);
266}
267
268MoFEMErrorCode Tools::writeTetsWithQuality(const char *file_name,
269 const char *file_type,
270 const char *options,
271 const Range &tets, Tag th,
272 boost::function<bool(double)> f) {
273 MoFEM::Interface &m_field = cOre;
274 moab::Interface &moab(m_field.get_moab());
275 MoFEMFunctionBegin;
276 Range out_tets;
277 CHKERR getTetsWithQuality(out_tets, tets, th, f);
278 EntityHandle meshset;
279 CHKERR moab.create_meshset(MESHSET_SET, meshset);
280 CHKERR moab.add_entities(meshset, out_tets);
281 CHKERR moab.write_file(file_name, file_type, options, &meshset, 1);
282 CHKERR moab.delete_entities(&meshset, 1);
283 MoFEMFunctionReturn(0);
284}
285
286MoFEMErrorCode Tools::checkIfPointIsInTet(const double tet_coords[],
287 const double global_coord[],
288 const double tol, bool &result) {
289 double loc_coord[] = {0, 0, 0};
290 double N[4], diffN[12];
291 MoFEMFunctionBegin;
292 CHKERR ShapeDiffMBTET(diffN);
293 CHKERR ShapeMBTET(N, &loc_coord[0], &loc_coord[1], &loc_coord[2], 1);
294 CHKERR ShapeMBTET_inverse(N, diffN, tet_coords, global_coord, loc_coord);
295 CHKERR ShapeMBTET(N, &loc_coord[0], &loc_coord[1], &loc_coord[2], 1);
296 result = true;
297 for (int n = 0; n != 4; ++n) {
298 if (N[n] < -tol || (N[n] - 1) > tol) {
299 result = false;
300 break;
301 }
302 }
303 MoFEMFunctionReturn(0);
304}
305
306MoFEMErrorCode Tools::checkVectorForNotANumber(const Problem *prb_ptr,
307 const RowColData row_or_col,
308 Vec v) {
309 MoFEMFunctionBegin;
310 int loc_size;
311 CHKERR VecGetLocalSize(v, &loc_size);
312 int prb_loc_size = 0;
313 boost::shared_ptr<NumeredDofEntity_multiIndex> prb_dofs;
314 switch (row_or_col) {
315 case ROW:
316 prb_loc_size = prb_ptr->getNbLocalDofsRow();
317 prb_dofs = prb_ptr->getNumeredRowDofsPtr();
318 break;
319 case COL:
320 prb_loc_size = prb_ptr->getNbLocalDofsCol();
321 prb_dofs = prb_ptr->getNumeredColDofsPtr();
322 break;
323 break;
324 default:
325 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
326 "Wrong argument, row_or_col should be row or column");
327 }
328 if (loc_size != prb_loc_size) {
329 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
330 "Inconsistent size of vector and problem %d != %d", loc_size,
331 prb_loc_size);
332 }
333 const double *a;
334 CHKERR VecGetArrayRead(v, &a);
335 MPI_Comm comm = PetscObjectComm((PetscObject)v);
336 for (int ii = 0; ii != loc_size; ++ii) {
337 if (!boost::math::isfinite(a[ii])) {
338 NumeredDofEntityByLocalIdx::iterator dit =
339 prb_dofs->get<PetscLocalIdx_mi_tag>().find(ii);
340 std::ostringstream ss;
341 ss << "Not a number " << a[ii] << " on dof: " << endl
342 << **dit << endl
343 << endl;
344 PetscSynchronizedPrintf(comm, "%s", ss.str().c_str());
345 }
346 }
347 CHKERR VecRestoreArrayRead(v, &a);
348 PetscSynchronizedFlush(comm, PETSC_STDOUT);
349 MoFEMFunctionReturn(0);
350}
351
352MoFEMErrorCode Tools::getTriNormal(const double *coords, double *normal) {
353 MoFEMFunctionBegin;
354 double diffN[6];
355 CHKERR ShapeDiffMBTRI(diffN);
356 CHKERR ShapeFaceNormalMBTRI(diffN, coords, normal);
357 MoFEMFunctionReturn(0);
358}
359
360MoFEMErrorCode Tools::getTriNormal(const EntityHandle tri,
361 double *normal) const {
362 MoFEM::Interface &m_field = cOre;
363 moab::Interface &moab(m_field.get_moab());
364 MoFEMFunctionBegin;
365 if (type_from_handle(tri) != MBTRI) {
366 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "Works only for triangle");
367 }
368 const EntityHandle *conn;
369 int num_nodes;
370 double coords[9];
371 CHKERR moab.get_connectivity(tri, conn, num_nodes, true);
372 CHKERR moab.get_coords(conn, num_nodes, coords);
373 CHKERR getTriNormal(coords, normal);
374 MoFEMFunctionReturn(0);
375}
376
377double Tools::getTriArea(const EntityHandle tri) const {
378 FTensor::Index<'i', 3> i;
379 FTensor::Tensor1<double, 3> t_normal;
380 CHK_THROW_MESSAGE(getTriNormal(tri, &t_normal(0)), "calculate area");
381 return sqrt(t_normal(i) * t_normal(i)) * 0.5;
382}
383
384double Tools::getEdgeLength(const double *edge_coords) {
385 FTensor::Tensor1<double, 3> t_coords_n0(edge_coords[0], edge_coords[1],
386 edge_coords[2]);
387 FTensor::Tensor1<double, 3> t_coords_n1(edge_coords[3], edge_coords[4],
388 edge_coords[5]);
389 FTensor::Index<'i', 3> i;
390 t_coords_n0(i) -= t_coords_n1(i);
391 return sqrt(t_coords_n0(i) * t_coords_n0(i));
392}
393
394double Tools::getEdgeLength(const EntityHandle edge) {
395 MoFEM::Interface &m_field = cOre;
396 moab::Interface &moab(m_field.get_moab());
397 auto get_edge_coords = [edge, &moab](double *const coords) {
398 MoFEMFunctionBegin;
399 if (type_from_handle(edge) != MBEDGE) {
400 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "Works only for edge");
401 }
402 const EntityHandle *conn;
403 int num_nodes;
404 CHKERR moab.get_connectivity(edge, conn, num_nodes, true);
405 CHKERR moab.get_coords(conn, 2, coords);
406 MoFEMFunctionReturn(0);
407 };
408 double coords[6];
409 ierr = get_edge_coords(coords);
410 CHKERRABORT(PETSC_COMM_SELF, ierr);
411 return getEdgeLength(coords);
412}
413
414Tools::SEGMENT_MIN_DISTANCE
415Tools::minDistancePointFromOnSegment(const double *w_ptr, const double *v_ptr,
416 const double *p_ptr, double *const t_ptr) {
417 FTensor::Index<'i', 3> i;
418 FTensor::Tensor1<const double *, 3> t_w(&w_ptr[0], &w_ptr[1], &w_ptr[2]);
419 FTensor::Tensor1<const double *, 3> t_v(&v_ptr[0], &v_ptr[1], &v_ptr[2]);
420 FTensor::Tensor1<const double *, 3> t_p(&p_ptr[0], &p_ptr[1], &p_ptr[2]);
421 FTensor::Tensor1<double, 3> t_vw;
422 t_vw(i) = t_v(i) - t_w(i);
423 const double dot_vw = t_vw(i) * t_vw(i);
424 if (std::abs(dot_vw) < std::numeric_limits<double>::epsilon()) {
425 if (t_ptr)
426 *t_ptr = 0;
427 return SEGMENT_ONE_IS_POINT;
428 }
429 FTensor::Tensor1<double, 3> t_pw;
430 t_pw(i) = t_p(i) - t_w(i);
431 const double t = t_pw(i) * t_vw(i) / dot_vw;
432 if (t_ptr)
433 *t_ptr = t;
434 return SOLUTION_EXIST;
435}
436
437Tools::SEGMENT_MIN_DISTANCE
438Tools::minDistanceFromSegments(const double *w_ptr, const double *v_ptr,
439 const double *k_ptr, const double *l_ptr,
440 double *const tvw_ptr, double *const tlk_ptr) {
441
442 FTensor::Index<'i', 3> i;
443 FTensor::Tensor1<const double *, 3> t_w(&w_ptr[0], &w_ptr[1], &w_ptr[2]);
444 FTensor::Tensor1<const double *, 3> t_v(&v_ptr[0], &v_ptr[1], &v_ptr[2]);
445 FTensor::Tensor1<const double *, 3> t_k(&k_ptr[0], &k_ptr[1], &k_ptr[2]);
446 FTensor::Tensor1<const double *, 3> t_l(&l_ptr[0], &l_ptr[1], &l_ptr[2]);
447
448 // First segnent is a point
449 FTensor::Tensor1<double, 3> t_vw;
450 t_vw(i) = t_v(i) - t_w(i);
451 double dot_vw = t_vw(i) * t_vw(i);
452 if (std::abs(dot_vw) < std::numeric_limits<double>::epsilon()) {
453 if (tvw_ptr)
454 *tvw_ptr = 0;
455 if (minDistancePointFromOnSegment(k_ptr, l_ptr, w_ptr, tlk_ptr) ==
456 SEGMENT_ONE_IS_POINT)
457 return SEGMENT_TWO_AND_TWO_ARE_POINT;
458 else
459 return SEGMENT_ONE_IS_POINT;
460 }
461
462 // Second segment is a point
463 FTensor::Tensor1<double, 3> t_lk;
464 t_lk(i) = t_l(i) - t_k(i);
465 double dot_lk = t_lk(i) * t_lk(i);
466 if (std::abs(dot_lk) < std::numeric_limits<double>::epsilon()) {
467 if (tlk_ptr)
468 *tlk_ptr = 0;
469 if (minDistancePointFromOnSegment(w_ptr, v_ptr, k_ptr, tvw_ptr) ==
470 SEGMENT_ONE_IS_POINT)
471 return SEGMENT_TWO_AND_TWO_ARE_POINT;
472 else
473 return SEGMENT_TWO_IS_POINT;
474 }
475
476 const double a = t_vw(i) * t_vw(i);
477 const double b = -t_vw(i) * t_lk(i);
478 const double c = t_lk(i) * t_lk(i);
479
480 const double det = a * c - b * b;
481 if (std::abs(det) < std::numeric_limits<double>::epsilon()) {
482
483 return NO_SOLUTION;
484
485 } else {
486
487 FTensor::Tensor1<double, 3> t_wk;
488 t_wk(i) = t_w(i) - t_k(i);
489
490 const double ft0 = t_vw(i) * t_wk(i);
491 const double ft1 = -t_lk(i) * t_wk(i);
492 const double t0 = (ft1 * b - ft0 * c) / det;
493 const double t1 = (ft0 * b - ft1 * a) / det;
494
495 if (tvw_ptr)
496 *tvw_ptr = t0;
497 if (tlk_ptr)
498 *tlk_ptr = t1;
499
500 return SOLUTION_EXIST;
501 }
502}
503
504MoFEMErrorCode Tools::findMinDistanceFromTheEdges(
505 const double *v_ptr, const int nb, Range edges, double *min_dist_ptr,
506 double *o_ptr, EntityHandle *o_segments) const {
507 MoFEM::Interface &m_field = cOre;
508 moab::Interface &moab(m_field.get_moab());
509 MoFEMFunctionBegin;
510
511 FTensor::Index<'i', 3> i;
512
513 auto get_point = [i](auto &t_w, auto &t_delta, auto t) {
514 FTensor::Tensor1<double, 3> t_p;
515 t = std::max(0., std::min(1., t));
516 t_p(i) = t_w(i) + t * t_delta(i);
517 return t_p;
518 };
519
520 auto get_distance = [i](auto &t_p, auto &t_n) {
521 FTensor::Tensor1<double, 3> t_dist_vector;
522 t_dist_vector(i) = t_p(i) - t_n(i);
523 return sqrt(t_dist_vector(i) * t_dist_vector(i));
524 };
525
526 for (auto e : edges) {
527
528 int num_nodes;
529 const EntityHandle *conn_fixed;
530 CHKERR moab.get_connectivity(e, conn_fixed, num_nodes, true);
531 VectorDouble6 coords_fixed(6);
532 CHKERR moab.get_coords(conn_fixed, num_nodes, &coords_fixed[0]);
533 FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_f0(
534 &coords_fixed[0], &coords_fixed[1], &coords_fixed[2]);
535 FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_f1(
536 &coords_fixed[3], &coords_fixed[4], &coords_fixed[5]);
537
538 FTensor::Tensor1<double, 3> t_edge_delta;
539 t_edge_delta(i) = t_f1(i) - t_f0(i);
540
541 FTensor::Tensor1<FTensor::PackPtr<const double *, 3>, 3> t_n(
542 v_ptr, v_ptr + 1, v_ptr + 2);
543 FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_min_coords(
544 o_ptr, o_ptr + 1, o_ptr + 2);
545 FTensor::Tensor0<FTensor::PackPtr<double *, 1>> t_min_dist(min_dist_ptr);
546
547 EntityHandle *colsest_segment_it = nullptr;
548 if (o_segments)
549 colsest_segment_it = o_segments;
550
551 for (int n = 0; n != nb; ++n) {
552
553 double t;
554 if (Tools::minDistancePointFromOnSegment(&t_f0(0), &t_f1(0), &t_n(0),
555 &t) == Tools::SOLUTION_EXIST) {
556 auto t_p = get_point(t_f0, t_edge_delta, t);
557 auto dist_n = get_distance(t_p, t_n);
558 if (dist_n < t_min_dist || t_min_dist < 0) {
559 t_min_dist = dist_n;
560 if (o_ptr)
561 t_min_coords(i) = t_p(i);
562 if (o_segments)
563 *colsest_segment_it = e;
564 }
565 }
566
567 ++t_n;
568 ++t_min_dist;
569 if (o_ptr)
570 ++t_min_coords;
571 if (o_segments)
572 ++colsest_segment_it;
573 }
574 }
575
576 MoFEMFunctionReturn(0);
577}
578
579MoFEMErrorCode Tools::outerProductOfEdgeIntegrationPtsForQuad(
580 MatrixDouble &gauss_pts, const int rule_ksi, const int rule_eta) {
581 MoFEMFunctionBegin;
582
583 auto check_rule_edge = [](int rule) {
584 MoFEMFunctionBeginHot;
585 if (rule < 0) {
586 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
587 "Wrong integration rule: %d", rule);
588 }
589 if (rule > QUAD_1D_TABLE_SIZE) {
590 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
591 "rule > quadrature order %d < %d", rule, QUAD_1D_TABLE_SIZE);
592 }
593 if (QUAD_1D_TABLE[rule]->dim != 1) {
594 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong dimension");
595 }
596 if (QUAD_1D_TABLE[rule]->order < rule) {
597 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
598 "wrong order %d != %d", QUAD_1D_TABLE[rule]->order, rule);
599 }
600 MoFEMFunctionReturnHot(0);
601 };
602
603 CHKERR check_rule_edge(rule_ksi);
604 CHKERR check_rule_edge(rule_eta);
605
606 const int nb_gauss_pts_ksi = QUAD_1D_TABLE[rule_ksi]->npoints;
607 const int nb_gauss_pts_eta = QUAD_1D_TABLE[rule_eta]->npoints;
608 gauss_pts.resize(3, nb_gauss_pts_ksi * nb_gauss_pts_eta, false);
609
610 int gg = 0;
611 for (size_t i = 0; i != nb_gauss_pts_ksi; ++i) {
612 const double wi = QUAD_1D_TABLE[rule_ksi]->weights[i];
613 const double ksi = (QUAD_1D_TABLE[rule_ksi]->points[2 * i + 1]);
614 for (size_t j = 0; j != nb_gauss_pts_eta; ++j, ++gg) {
615 const double wk = wi * QUAD_1D_TABLE[rule_eta]->weights[j];
616 const double eta = QUAD_1D_TABLE[rule_eta]->points[2 * j + 1];
617 gauss_pts(0, gg) = ksi;
618 gauss_pts(1, gg) = eta;
619 gauss_pts(2, gg) = wk;
620 }
621 }
622 if (gg != gauss_pts.size2())
623 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong size of matrix");
624
625 MoFEMFunctionReturn(0);
626}
627
628MoFEMErrorCode Tools::outerProductOfEdgeIntegrationPtsForHex(
629 MatrixDouble &gauss_pts, const int rule_ksi, const int rule_eta,
630 const int rule_zeta) {
631 MoFEMFunctionBegin;
632
633 auto check_rule_edge = [](int rule) {
634 MoFEMFunctionBeginHot;
635 if (rule < 0) {
636 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
637 "Wrong integration rule: %d", rule);
638 }
639 if (rule > QUAD_1D_TABLE_SIZE) {
640 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
641 "rule > quadrature order %d < %d", rule, QUAD_1D_TABLE_SIZE);
642 }
643 if (QUAD_1D_TABLE[rule]->dim != 1) {
644 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong dimension");
645 }
646 if (QUAD_1D_TABLE[rule]->order < rule) {
647 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
648 "wrong order %d != %d", QUAD_1D_TABLE[rule]->order, rule);
649 }
650 MoFEMFunctionReturnHot(0);
651 };
652
653 CHKERR check_rule_edge(rule_ksi);
654 CHKERR check_rule_edge(rule_eta);
655 CHKERR check_rule_edge(rule_zeta);
656
657 const int nb_gauss_pts_ksi = QUAD_1D_TABLE[rule_ksi]->npoints;
658 const int nb_gauss_pts_eta = QUAD_1D_TABLE[rule_eta]->npoints;
659 const int nb_gauss_pts_zeta = QUAD_1D_TABLE[rule_zeta]->npoints;
660 gauss_pts.resize(4, nb_gauss_pts_ksi * nb_gauss_pts_eta * nb_gauss_pts_zeta,
661 false);
662
663 int gg = 0;
664 for (size_t i = 0; i != nb_gauss_pts_ksi; ++i) {
665 const double wi = QUAD_1D_TABLE[rule_ksi]->weights[i];
666 const double ksi = QUAD_1D_TABLE[rule_ksi]->points[2 * i + 1];
667 for (size_t j = 0; j != nb_gauss_pts_eta; ++j) {
668 const double wj = wi * QUAD_1D_TABLE[rule_eta]->weights[j];
669 const double eta = QUAD_1D_TABLE[rule_eta]->points[2 * j + 1];
670 for (size_t k = 0; k != nb_gauss_pts_zeta; ++k, ++gg) {
671 const double wk = wj * QUAD_1D_TABLE[rule_zeta]->weights[k];
672 const double zeta = QUAD_1D_TABLE[rule_zeta]->points[2 * k + 1];
673 gauss_pts(0, gg) = ksi;
674 gauss_pts(1, gg) = eta;
675 gauss_pts(2, gg) = zeta;
676 gauss_pts(3, gg) = wk;
677 }
678 }
679 }
680
681 if (gg != gauss_pts.size2())
682 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong size of matrix");
683
684 MoFEMFunctionReturn(0);
685}
686
687} // namespace MoFEM
a
constexpr double a
Definition: approx_sphere.cpp:30
FTensor::Tensor2_symmetric
Definition: Tensor2_symmetric_value.hpp:14
RowColData
RowColData
RowColData.
Definition: definitions.h:123
COL
@ COL
Definition: definitions.h:123
ROW
@ ROW
Definition: definitions.h:123
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:595
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
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:346
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:416
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
dim
const int dim
Definition: elec_phys_2D.cpp:12
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
ShapeFaceNormalMBTRI
PetscErrorCode ShapeFaceNormalMBTRI(double *diffN, const double *coords, double *normal)
Definition: fem_tools.c:229
ShapeDiffMBTRI
PetscErrorCode ShapeDiffMBTRI(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:194
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
n
FTensor::Index< 'n', SPACE_DIM > n
Definition: fluid_structure_eigenproblem.cpp:66
order
int order
Definition: fluid_structure_eigenproblem.cpp:84
edge_coords
static const double edge_coords[6][6]
Definition: forces_and_sources_h1_continuity_check.cpp:18
eta
constexpr double eta
Definition: free_surface.cpp:96
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
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
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:27
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
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:1634
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1352
zeta
double zeta
Definition: plastic.cpp:107
t
constexpr double t
plate stiffness
Definition: plate.cpp:59
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
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:40
MoFEM::Problem::getNbLocalDofsRow
DofIdx getNbLocalDofsRow() const
Definition: ProblemsMultiIndices.hpp:356
MoFEM::Problem::getNbLocalDofsCol
DofIdx getNbLocalDofsCol() const
Definition: ProblemsMultiIndices.hpp:357
MoFEM::Problem::getNumeredColDofsPtr
auto & getNumeredColDofsPtr() const
get access to numeredColDofsPtr storing DOFs on cols
Definition: ProblemsMultiIndices.hpp:73
MoFEM::Problem::getNumeredRowDofsPtr
auto & getNumeredRowDofsPtr() const
get access to numeredRowDofsPtr storing DOFs on rows
Definition: ProblemsMultiIndices.hpp:68
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:438
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:90
MoFEM::Tools::outerProductOfEdgeIntegrationPtsForQuad
static MoFEMErrorCode outerProductOfEdgeIntegrationPtsForQuad(MatrixDouble &pts, const int edge0, const int edge1)
Definition: Tools.cpp:579
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:415
MoFEM::Tools::shapeFunMBTRI1At00
static constexpr double shapeFunMBTRI1At00
Definition: Tools.hpp:113
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:504
MoFEM::Tools::diffShapeFunMBTRI
static constexpr std::array< double, 6 > diffShapeFunMBTRI
Definition: Tools.hpp:104
MoFEM::Tools::SEGMENT_MIN_DISTANCE
SEGMENT_MIN_DISTANCE
Definition: Tools.hpp:507
MoFEM::Tools::SOLUTION_EXIST
@ SOLUTION_EXIST
Definition: Tools.hpp:508
MoFEM::Tools::SEGMENT_TWO_AND_TWO_ARE_POINT
@ SEGMENT_TWO_AND_TWO_ARE_POINT
Definition: Tools.hpp:511
MoFEM::Tools::SEGMENT_ONE_IS_POINT
@ SEGMENT_ONE_IS_POINT
Definition: Tools.hpp:509
MoFEM::Tools::SEGMENT_TWO_IS_POINT
@ SEGMENT_TWO_IS_POINT
Definition: Tools.hpp:510
MoFEM::Tools::getLocalCoordinatesOnReferenceTriNodeTri
static MoFEMErrorCode getLocalCoordinatesOnReferenceTriNodeTri(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:142
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:306
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::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:51
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:268
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:243
MoFEM::Tools::outerProductOfEdgeIntegrationPtsForHex
static MoFEMErrorCode outerProductOfEdgeIntegrationPtsForHex(MatrixDouble &pts, const int edge0, const int edge1, const int edge2)
Definition: Tools.cpp:628
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:286
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:200
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:32
MoFEM::Tools::volumeLengthQuality
static double volumeLengthQuality(const double *coords)
Calculate tetrahedron volume length quality.
Definition: Tools.cpp:17
MoFEM::Tools::getTriArea
double getTriArea(const EntityHandle tri) const
Get triangle area.
Definition: Tools.cpp:377
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:11
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:384
MoFEM::Tools::cOre
MoFEM::Core & cOre
Definition: Tools.hpp:24
MoFEM::Tools::getTriNormal
static MoFEMErrorCode getTriNormal(const double *coords, double *normal)
Get the Tri Normal objectGet triangle normal.
Definition: Tools.cpp:352
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.9.5 and hosted at University of Glasgow