v0.14.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
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<double, 3> t_coords_at_0;
211 // Build matrix and get coordinates of zero point
212 // ii - global coordinates
213 FTensor::Tensor1<double, 3> t_a;
214 for (auto ii : {0, 1, 2}) {
215 t_a(ii) = diffShapeFunMBEDGE[0] * t_elem_coords(0) +
216 diffShapeFunMBEDGE[1] * t_elem_coords(1);
217 t_coords_at_0(ii) = shapeFunMBEDGEAt00[0] * t_elem_coords(0) +
218 shapeFunMBEDGEAt00[1] * t_elem_coords(1);
219 ++t_elem_coords;
220 }
221
222 FTensor::Tensor1<FTensor::PackPtr<const double *, 3>, 3> t_global_coords = {
223 &global_coords[0], &global_coords[1], &global_coords[2]};
224 FTensor::Tensor0<FTensor::PackPtr<double *, 1>> t_local_coords = {
225 &local_coords[0]};
226
227 const double b = t_a(i3) * t_a(i3);
228 const double inv_b = 1 / b;
229
230 // Calculate right hand side
231 for (int ii = 0; ii != nb_nodes; ++ii) {
232 t_local_coords =
233 inv_b * (t_a(i3) * (t_global_coords(i3) - t_coords_at_0(i3)));
234 ++t_local_coords;
235 ++t_global_coords;
236 }
237
238 MoFEMFunctionReturnHot(0);
239}
240
241MoFEMErrorCode Tools::getTetsWithQuality(Range &out_tets, const Range &tets,
242 Tag th,
243 boost::function<bool(double)> f) {
244 MoFEM::Interface &m_field = cOre;
245 moab::Interface &moab(m_field.get_moab());
246 MoFEMFunctionBegin;
247 Range to_write;
248 const EntityHandle *conn;
249 int num_nodes;
250 double coords[12];
251 for (auto tet : tets) {
252 CHKERR m_field.get_moab().get_connectivity(tet, conn, num_nodes, true);
253 if (th) {
254 CHKERR moab.tag_get_data(th, conn, num_nodes, coords);
255 } else {
256 CHKERR moab.get_coords(conn, num_nodes, coords);
257 }
258 double q = Tools::volumeLengthQuality(coords);
259 if (f(q)) {
260 out_tets.insert(tet);
261 }
262 }
263 MoFEMFunctionReturn(0);
264}
265
266MoFEMErrorCode Tools::writeTetsWithQuality(const char *file_name,
267 const char *file_type,
268 const char *options,
269 const Range &tets, Tag th,
270 boost::function<bool(double)> f) {
271 MoFEM::Interface &m_field = cOre;
272 moab::Interface &moab(m_field.get_moab());
273 MoFEMFunctionBegin;
274 Range out_tets;
275 CHKERR getTetsWithQuality(out_tets, tets, th, f);
276 EntityHandle meshset;
277 CHKERR moab.create_meshset(MESHSET_SET, meshset);
278 CHKERR moab.add_entities(meshset, out_tets);
279 CHKERR moab.write_file(file_name, file_type, options, &meshset, 1);
280 CHKERR moab.delete_entities(&meshset, 1);
281 MoFEMFunctionReturn(0);
282}
283
284MoFEMErrorCode Tools::checkIfPointIsInTet(const double tet_coords[],
285 const double global_coord[],
286 const double tol, bool &result) {
287 double loc_coord[] = {0, 0, 0};
288 double N[4], diffN[12];
289 MoFEMFunctionBegin;
290 CHKERR ShapeDiffMBTET(diffN);
291 CHKERR ShapeMBTET(N, &loc_coord[0], &loc_coord[1], &loc_coord[2], 1);
292 CHKERR ShapeMBTET_inverse(N, diffN, tet_coords, global_coord, loc_coord);
293 CHKERR ShapeMBTET(N, &loc_coord[0], &loc_coord[1], &loc_coord[2], 1);
294 result = true;
295 for (int n = 0; n != 4; ++n) {
296 if (N[n] < -tol || (N[n] - 1) > tol) {
297 result = false;
298 break;
299 }
300 }
301 MoFEMFunctionReturn(0);
302}
303
304MoFEMErrorCode Tools::checkVectorForNotANumber(const Problem *prb_ptr,
305 const RowColData row_or_col,
306 Vec v) {
307 MoFEMFunctionBegin;
308 int loc_size;
309 CHKERR VecGetLocalSize(v, &loc_size);
310 int prb_loc_size = 0;
311 boost::shared_ptr<NumeredDofEntity_multiIndex> prb_dofs;
312 switch (row_or_col) {
313 case ROW:
314 prb_loc_size = prb_ptr->getNbLocalDofsRow();
315 prb_dofs = prb_ptr->getNumeredRowDofsPtr();
316 break;
317 case COL:
318 prb_loc_size = prb_ptr->getNbLocalDofsCol();
319 prb_dofs = prb_ptr->getNumeredColDofsPtr();
320 break;
321 break;
322 default:
323 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
324 "Wrong argument, row_or_col should be row or column");
325 }
326 if (loc_size != prb_loc_size) {
327 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
328 "Inconsistent size of vector and problem %d != %d", loc_size,
329 prb_loc_size);
330 }
331 const double *a;
332 CHKERR VecGetArrayRead(v, &a);
333 MPI_Comm comm = PetscObjectComm((PetscObject)v);
334 for (int ii = 0; ii != loc_size; ++ii) {
335 if (!boost::math::isfinite(a[ii])) {
336 NumeredDofEntityByLocalIdx::iterator dit =
337 prb_dofs->get<PetscLocalIdx_mi_tag>().find(ii);
338 std::ostringstream ss;
339 ss << "Not a number " << a[ii] << " on dof: " << endl
340 << **dit << endl
341 << endl;
342 PetscSynchronizedPrintf(comm, "%s", ss.str().c_str());
343 }
344 }
345 CHKERR VecRestoreArrayRead(v, &a);
346 PetscSynchronizedFlush(comm, PETSC_STDOUT);
347 MoFEMFunctionReturn(0);
348}
349
350MoFEMErrorCode Tools::getTriNormal(const double *coords, double *normal) {
351 MoFEMFunctionBegin;
352 double diffN[6];
353 CHKERR ShapeDiffMBTRI(diffN);
354 CHKERR ShapeFaceNormalMBTRI(diffN, coords, normal);
355 MoFEMFunctionReturn(0);
356}
357
358MoFEMErrorCode Tools::getTriNormal(const EntityHandle tri,
359 double *normal) const {
360 MoFEM::Interface &m_field = cOre;
361 moab::Interface &moab(m_field.get_moab());
362 MoFEMFunctionBegin;
363 if (type_from_handle(tri) != MBTRI) {
364 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "Works only for triangle");
365 }
366 const EntityHandle *conn;
367 int num_nodes;
368 double coords[9];
369 CHKERR moab.get_connectivity(tri, conn, num_nodes, true);
370 CHKERR moab.get_coords(conn, num_nodes, coords);
371 CHKERR getTriNormal(coords, normal);
372 MoFEMFunctionReturn(0);
373}
374
375double Tools::getTriArea(const EntityHandle tri) const {
376 FTensor::Index<'i', 3> i;
377 FTensor::Tensor1<double, 3> t_normal;
378 CHK_THROW_MESSAGE(getTriNormal(tri, &t_normal(0)), "calculate area");
379 return sqrt(t_normal(i) * t_normal(i)) * 0.5;
380}
381
382double Tools::getEdgeLength(const double *edge_coords) {
383 FTensor::Tensor1<double, 3> t_coords_n0(edge_coords[0], edge_coords[1],
384 edge_coords[2]);
385 FTensor::Tensor1<double, 3> t_coords_n1(edge_coords[3], edge_coords[4],
386 edge_coords[5]);
387 FTensor::Index<'i', 3> i;
388 t_coords_n0(i) -= t_coords_n1(i);
389 return sqrt(t_coords_n0(i) * t_coords_n0(i));
390}
391
392double Tools::getEdgeLength(const EntityHandle edge) {
393 MoFEM::Interface &m_field = cOre;
394 moab::Interface &moab(m_field.get_moab());
395 auto get_edge_coords = [edge, &moab](double *const coords) {
396 MoFEMFunctionBegin;
397 if (type_from_handle(edge) != MBEDGE) {
398 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "Works only for edge");
399 }
400 const EntityHandle *conn;
401 int num_nodes;
402 CHKERR moab.get_connectivity(edge, conn, num_nodes, true);
403 CHKERR moab.get_coords(conn, 2, coords);
404 MoFEMFunctionReturn(0);
405 };
406 double coords[6];
407 ierr = get_edge_coords(coords);
408 CHKERRABORT(PETSC_COMM_SELF, ierr);
409 return getEdgeLength(coords);
410}
411
412Tools::SEGMENT_MIN_DISTANCE
413Tools::minDistancePointFromOnSegment(const double *w_ptr, const double *v_ptr,
414 const double *p_ptr, double *const t_ptr) {
415 FTensor::Index<'i', 3> i;
416 FTensor::Tensor1<const double *, 3> t_w(&w_ptr[0], &w_ptr[1], &w_ptr[2]);
417 FTensor::Tensor1<const double *, 3> t_v(&v_ptr[0], &v_ptr[1], &v_ptr[2]);
418 FTensor::Tensor1<const double *, 3> t_p(&p_ptr[0], &p_ptr[1], &p_ptr[2]);
419 FTensor::Tensor1<double, 3> t_vw;
420 t_vw(i) = t_v(i) - t_w(i);
421 const double dot_vw = t_vw(i) * t_vw(i);
422 if (std::abs(dot_vw) < std::numeric_limits<double>::epsilon()) {
423 if (t_ptr)
424 *t_ptr = 0;
425 return SEGMENT_ONE_IS_POINT;
426 }
427 FTensor::Tensor1<double, 3> t_pw;
428 t_pw(i) = t_p(i) - t_w(i);
429 const double t = t_pw(i) * t_vw(i) / dot_vw;
430 if (t_ptr)
431 *t_ptr = t;
432 return SOLUTION_EXIST;
433}
434
435Tools::SEGMENT_MIN_DISTANCE
436Tools::minDistanceFromSegments(const double *w_ptr, const double *v_ptr,
437 const double *k_ptr, const double *l_ptr,
438 double *const tvw_ptr, double *const tlk_ptr) {
439
440 FTensor::Index<'i', 3> i;
441 FTensor::Tensor1<const double *, 3> t_w(&w_ptr[0], &w_ptr[1], &w_ptr[2]);
442 FTensor::Tensor1<const double *, 3> t_v(&v_ptr[0], &v_ptr[1], &v_ptr[2]);
443 FTensor::Tensor1<const double *, 3> t_k(&k_ptr[0], &k_ptr[1], &k_ptr[2]);
444 FTensor::Tensor1<const double *, 3> t_l(&l_ptr[0], &l_ptr[1], &l_ptr[2]);
445
446 // First segnent is a point
447 FTensor::Tensor1<double, 3> t_vw;
448 t_vw(i) = t_v(i) - t_w(i);
449 double dot_vw = t_vw(i) * t_vw(i);
450 if (std::abs(dot_vw) < std::numeric_limits<double>::epsilon()) {
451 if (tvw_ptr)
452 *tvw_ptr = 0;
453 if (minDistancePointFromOnSegment(k_ptr, l_ptr, w_ptr, tlk_ptr) ==
454 SEGMENT_ONE_IS_POINT)
455 return SEGMENT_TWO_AND_TWO_ARE_POINT;
456 else
457 return SEGMENT_ONE_IS_POINT;
458 }
459
460 // Second segment is a point
461 FTensor::Tensor1<double, 3> t_lk;
462 t_lk(i) = t_l(i) - t_k(i);
463 double dot_lk = t_lk(i) * t_lk(i);
464 if (std::abs(dot_lk) < std::numeric_limits<double>::epsilon()) {
465 if (tlk_ptr)
466 *tlk_ptr = 0;
467 if (minDistancePointFromOnSegment(w_ptr, v_ptr, k_ptr, tvw_ptr) ==
468 SEGMENT_ONE_IS_POINT)
469 return SEGMENT_TWO_AND_TWO_ARE_POINT;
470 else
471 return SEGMENT_TWO_IS_POINT;
472 }
473
474 const double a = t_vw(i) * t_vw(i);
475 const double b = -t_vw(i) * t_lk(i);
476 const double c = t_lk(i) * t_lk(i);
477
478 const double det = a * c - b * b;
479 if (std::abs(det) < std::numeric_limits<double>::epsilon()) {
480
481 return NO_SOLUTION;
482
483 } else {
484
485 FTensor::Tensor1<double, 3> t_wk;
486 t_wk(i) = t_w(i) - t_k(i);
487
488 const double ft0 = t_vw(i) * t_wk(i);
489 const double ft1 = -t_lk(i) * t_wk(i);
490 const double t0 = (ft1 * b - ft0 * c) / det;
491 const double t1 = (ft0 * b - ft1 * a) / det;
492
493 if (tvw_ptr)
494 *tvw_ptr = t0;
495 if (tlk_ptr)
496 *tlk_ptr = t1;
497
498 return SOLUTION_EXIST;
499 }
500}
501
502MoFEMErrorCode Tools::findMinDistanceFromTheEdges(
503 const double *v_ptr, const int nb, Range edges, double *min_dist_ptr,
504 double *o_ptr, EntityHandle *o_segments) const {
505 MoFEM::Interface &m_field = cOre;
506 moab::Interface &moab(m_field.get_moab());
507 MoFEMFunctionBegin;
508
509 FTensor::Index<'i', 3> i;
510
511 auto get_point = [i](auto &t_w, auto &t_delta, auto t) {
512 FTensor::Tensor1<double, 3> t_p;
513 t = std::max(0., std::min(1., t));
514 t_p(i) = t_w(i) + t * t_delta(i);
515 return t_p;
516 };
517
518 auto get_distance = [i](auto &t_p, auto &t_n) {
519 FTensor::Tensor1<double, 3> t_dist_vector;
520 t_dist_vector(i) = t_p(i) - t_n(i);
521 return sqrt(t_dist_vector(i) * t_dist_vector(i));
522 };
523
524 for (auto e : edges) {
525
526 int num_nodes;
527 const EntityHandle *conn_fixed;
528 CHKERR moab.get_connectivity(e, conn_fixed, num_nodes, true);
529 VectorDouble6 coords_fixed(6);
530 CHKERR moab.get_coords(conn_fixed, num_nodes, &coords_fixed[0]);
531 FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_f0(
532 &coords_fixed[0], &coords_fixed[1], &coords_fixed[2]);
533 FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_f1(
534 &coords_fixed[3], &coords_fixed[4], &coords_fixed[5]);
535
536 FTensor::Tensor1<double, 3> t_edge_delta;
537 t_edge_delta(i) = t_f1(i) - t_f0(i);
538
539 FTensor::Tensor1<FTensor::PackPtr<const double *, 3>, 3> t_n(
540 v_ptr, v_ptr + 1, v_ptr + 2);
541 FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_min_coords(
542 o_ptr, o_ptr + 1, o_ptr + 2);
543 FTensor::Tensor0<FTensor::PackPtr<double *, 1>> t_min_dist(min_dist_ptr);
544
545 EntityHandle *colsest_segment_it = nullptr;
546 if (o_segments)
547 colsest_segment_it = o_segments;
548
549 for (int n = 0; n != nb; ++n) {
550
551 double t;
552 if (Tools::minDistancePointFromOnSegment(&t_f0(0), &t_f1(0), &t_n(0),
553 &t) == Tools::SOLUTION_EXIST) {
554 auto t_p = get_point(t_f0, t_edge_delta, t);
555 auto dist_n = get_distance(t_p, t_n);
556 if (dist_n < t_min_dist || t_min_dist < 0) {
557 t_min_dist = dist_n;
558 if (o_ptr)
559 t_min_coords(i) = t_p(i);
560 if (o_segments)
561 *colsest_segment_it = e;
562 }
563 }
564
565 ++t_n;
566 ++t_min_dist;
567 if (o_ptr)
568 ++t_min_coords;
569 if (o_segments)
570 ++colsest_segment_it;
571 }
572 }
573
574 MoFEMFunctionReturn(0);
575}
576
577MoFEMErrorCode Tools::outerProductOfEdgeIntegrationPtsForQuad(
578 MatrixDouble &gauss_pts, const int rule_ksi, const int rule_eta) {
579 MoFEMFunctionBegin;
580
581 auto check_rule_edge = [](int rule) {
582 MoFEMFunctionBeginHot;
583 if (rule < 0) {
584 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
585 "Wrong integration rule: %d", rule);
586 }
587 if (rule > QUAD_1D_TABLE_SIZE) {
588 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
589 "rule > quadrature order %d < %d", rule, QUAD_1D_TABLE_SIZE);
590 }
591 if (QUAD_1D_TABLE[rule]->dim != 1) {
592 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong dimension");
593 }
594 if (QUAD_1D_TABLE[rule]->order < rule) {
595 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
596 "wrong order %d != %d", QUAD_1D_TABLE[rule]->order, rule);
597 }
598 MoFEMFunctionReturnHot(0);
599 };
600
601 CHKERR check_rule_edge(rule_ksi);
602 CHKERR check_rule_edge(rule_eta);
603
604 const int nb_gauss_pts_ksi = QUAD_1D_TABLE[rule_ksi]->npoints;
605 const int nb_gauss_pts_eta = QUAD_1D_TABLE[rule_eta]->npoints;
606 gauss_pts.resize(3, nb_gauss_pts_ksi * nb_gauss_pts_eta, false);
607
608 int gg = 0;
609 for (size_t i = 0; i != nb_gauss_pts_ksi; ++i) {
610 const double wi = QUAD_1D_TABLE[rule_ksi]->weights[i];
611 const double ksi = (QUAD_1D_TABLE[rule_ksi]->points[2 * i + 1]);
612 for (size_t j = 0; j != nb_gauss_pts_eta; ++j, ++gg) {
613 const double wk = wi * QUAD_1D_TABLE[rule_eta]->weights[j];
614 const double eta = QUAD_1D_TABLE[rule_eta]->points[2 * j + 1];
615 gauss_pts(0, gg) = ksi;
616 gauss_pts(1, gg) = eta;
617 gauss_pts(2, gg) = wk;
618 }
619 }
620 if (gg != gauss_pts.size2())
621 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong size of matrix");
622
623 MoFEMFunctionReturn(0);
624}
625
626MoFEMErrorCode Tools::outerProductOfEdgeIntegrationPtsForHex(
627 MatrixDouble &gauss_pts, const int rule_ksi, const int rule_eta,
628 const int rule_zeta) {
629 MoFEMFunctionBegin;
630
631 auto check_rule_edge = [](int rule) {
632 MoFEMFunctionBeginHot;
633 if (rule < 0) {
634 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
635 "Wrong integration rule: %d", rule);
636 }
637 if (rule > QUAD_1D_TABLE_SIZE) {
638 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
639 "rule > quadrature order %d < %d", rule, QUAD_1D_TABLE_SIZE);
640 }
641 if (QUAD_1D_TABLE[rule]->dim != 1) {
642 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong dimension");
643 }
644 if (QUAD_1D_TABLE[rule]->order < rule) {
645 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
646 "wrong order %d != %d", QUAD_1D_TABLE[rule]->order, rule);
647 }
648 MoFEMFunctionReturnHot(0);
649 };
650
651 CHKERR check_rule_edge(rule_ksi);
652 CHKERR check_rule_edge(rule_eta);
653 CHKERR check_rule_edge(rule_zeta);
654
655 const int nb_gauss_pts_ksi = QUAD_1D_TABLE[rule_ksi]->npoints;
656 const int nb_gauss_pts_eta = QUAD_1D_TABLE[rule_eta]->npoints;
657 const int nb_gauss_pts_zeta = QUAD_1D_TABLE[rule_zeta]->npoints;
658 gauss_pts.resize(4, nb_gauss_pts_ksi * nb_gauss_pts_eta * nb_gauss_pts_zeta,
659 false);
660
661 int gg = 0;
662 for (size_t i = 0; i != nb_gauss_pts_ksi; ++i) {
663 const double wi = QUAD_1D_TABLE[rule_ksi]->weights[i];
664 const double ksi = QUAD_1D_TABLE[rule_ksi]->points[2 * i + 1];
665 for (size_t j = 0; j != nb_gauss_pts_eta; ++j) {
666 const double wj = wi * QUAD_1D_TABLE[rule_eta]->weights[j];
667 const double eta = QUAD_1D_TABLE[rule_eta]->points[2 * j + 1];
668 for (size_t k = 0; k != nb_gauss_pts_zeta; ++k, ++gg) {
669 const double wk = wj * QUAD_1D_TABLE[rule_zeta]->weights[k];
670 const double zeta = QUAD_1D_TABLE[rule_zeta]->points[2 * k + 1];
671 gauss_pts(0, gg) = ksi;
672 gauss_pts(1, gg) = eta;
673 gauss_pts(2, gg) = zeta;
674 gauss_pts(3, gg) = wk;
675 }
676 }
677 }
678
679 if (gg != gauss_pts.size2())
680 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong size of matrix");
681
682 MoFEMFunctionReturn(0);
683}
684
685} // 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
double eta
Definition: free_surface.cpp:170
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:1761
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1475
zeta
double zeta
Viscous hardening.
Definition: plastic.cpp:170
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: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:436
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:577
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:413
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:502
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:304
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:266
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:241
MoFEM::Tools::outerProductOfEdgeIntegrationPtsForHex
static MoFEMErrorCode outerProductOfEdgeIntegrationPtsForHex(MatrixDouble &pts, const int edge0, const int edge1, const int edge2)
Definition: Tools.cpp:626
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:284
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:375
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:382
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:350
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