v0.13.1
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);
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
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());
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 }
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
91 const double *elem_coords, const double *global_coords, const int nb_nodes,
92 double *local_coords) {
95
97 &elem_coords[0], &elem_coords[3], &elem_coords[6], &elem_coords[9]};
98
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}) {
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
121 &global_coords[0], &global_coords[1], &global_coords[2]};
123 &local_coords[0], &local_coords[1], &local_coords[2]};
124
125 // Calculate right hand side
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
140}
141
143 const double *elem_coords, const double *global_coords, const int nb_nodes,
144 double *local_coords) {
145
151
153 &elem_coords[0], &elem_coords[3], &elem_coords[6]};
154
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
163 for (auto ii : {0, 1, 2}) {
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
175 &global_coords[0], &global_coords[1], &global_coords[2]};
177 &local_coords[0], &local_coords[1]};
178
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
198}
199
201 const double *elem_coords, const double *global_coords, const int nb_nodes,
202 double *local_coords) {
203
206
208 &elem_coords[0], &elem_coords[3], &elem_coords[6]};
209
212 FTensor::Tensor1<double, 3> t_coords_at_0;
213 // Build matrix and get coordinates of zero point
214 // ii - global coordinates
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
225 &global_coords[0], &global_coords[1], &global_coords[2]};
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
241}
242
244 Tag th,
245 boost::function<bool(double)> f) {
246 MoFEM::Interface &m_field = cOre;
247 moab::Interface &moab(m_field.get_moab());
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 }
266}
267
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());
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);
284}
285
287 const double global_coord[],
288 const double tol, bool &result) {
289 double loc_coord[] = {0, 0, 0};
290 double N[4], diffN[12];
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 }
304}
305
307 const RowColData row_or_col,
308 Vec v) {
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);
350}
351
352MoFEMErrorCode Tools::getTriNormal(const double *coords, double *normal) {
354 double diffN[6];
355 CHKERR ShapeDiffMBTRI(diffN);
356 CHKERR ShapeFaceNormalMBTRI(diffN, coords, normal);
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());
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);
375}
376
377double Tools::getTriArea(const EntityHandle tri) const {
379 ierr = getTriNormal(tri, &t_normal(0));
380 CHKERRABORT(PETSC_COMM_SELF, ierr);
382 return sqrt(t_normal(i) * t_normal(i)) * 0.5;
383}
384
385double Tools::getEdgeLength(const double *edge_coords) {
387 edge_coords[2]);
389 edge_coords[5]);
391 t_coords_n0(i) -= t_coords_n1(i);
392 return sqrt(t_coords_n0(i) * t_coords_n0(i));
393}
394
395double Tools::getEdgeLength(const EntityHandle edge) {
396 MoFEM::Interface &m_field = cOre;
397 moab::Interface &moab(m_field.get_moab());
398 auto get_edge_coords = [edge, &moab](double *const coords) {
400 if (type_from_handle(edge) != MBEDGE) {
401 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "Works only for edge");
402 }
403 const EntityHandle *conn;
404 int num_nodes;
405 CHKERR moab.get_connectivity(edge, conn, num_nodes, true);
406 CHKERR moab.get_coords(conn, 2, coords);
408 };
409 double coords[6];
410 ierr = get_edge_coords(coords);
411 CHKERRABORT(PETSC_COMM_SELF, ierr);
412 return getEdgeLength(coords);
413}
414
416Tools::minDistancePointFromOnSegment(const double *w_ptr, const double *v_ptr,
417 const double *p_ptr, double *const t_ptr) {
419 FTensor::Tensor1<const double *, 3> t_w(&w_ptr[0], &w_ptr[1], &w_ptr[2]);
420 FTensor::Tensor1<const double *, 3> t_v(&v_ptr[0], &v_ptr[1], &v_ptr[2]);
421 FTensor::Tensor1<const double *, 3> t_p(&p_ptr[0], &p_ptr[1], &p_ptr[2]);
423 t_vw(i) = t_v(i) - t_w(i);
424 const double dot_vw = t_vw(i) * t_vw(i);
425 if (std::abs(dot_vw) < std::numeric_limits<double>::epsilon()) {
426 if (t_ptr)
427 *t_ptr = 0;
429 }
431 t_pw(i) = t_p(i) - t_w(i);
432 const double t = t_pw(i) * t_vw(i) / dot_vw;
433 if (t_ptr)
434 *t_ptr = t;
435 return SOLUTION_EXIST;
436}
437
439Tools::minDistanceFromSegments(const double *w_ptr, const double *v_ptr,
440 const double *k_ptr, const double *l_ptr,
441 double *const tvw_ptr, double *const tlk_ptr) {
442
444 FTensor::Tensor1<const double *, 3> t_w(&w_ptr[0], &w_ptr[1], &w_ptr[2]);
445 FTensor::Tensor1<const double *, 3> t_v(&v_ptr[0], &v_ptr[1], &v_ptr[2]);
446 FTensor::Tensor1<const double *, 3> t_k(&k_ptr[0], &k_ptr[1], &k_ptr[2]);
447 FTensor::Tensor1<const double *, 3> t_l(&l_ptr[0], &l_ptr[1], &l_ptr[2]);
448
449 // First segnent is a point
451 t_vw(i) = t_v(i) - t_w(i);
452 double dot_vw = t_vw(i) * t_vw(i);
453 if (std::abs(dot_vw) < std::numeric_limits<double>::epsilon()) {
454 if (tvw_ptr)
455 *tvw_ptr = 0;
456 if (minDistancePointFromOnSegment(k_ptr, l_ptr, w_ptr, tlk_ptr) ==
459 else
461 }
462
463 // Second segment is a point
465 t_lk(i) = t_l(i) - t_k(i);
466 double dot_lk = t_lk(i) * t_lk(i);
467 if (std::abs(dot_lk) < std::numeric_limits<double>::epsilon()) {
468 if (tlk_ptr)
469 *tlk_ptr = 0;
470 if (minDistancePointFromOnSegment(w_ptr, v_ptr, k_ptr, tvw_ptr) ==
473 else
475 }
476
477 const double a = t_vw(i) * t_vw(i);
478 const double b = -t_vw(i) * t_lk(i);
479 const double c = t_lk(i) * t_lk(i);
480
481 const double det = a * c - b * b;
482 if (std::abs(det) < std::numeric_limits<double>::epsilon()) {
483
484 return NO_SOLUTION;
485
486 } else {
487
489 t_wk(i) = t_w(i) - t_k(i);
490
491 const double ft0 = t_vw(i) * t_wk(i);
492 const double ft1 = -t_lk(i) * t_wk(i);
493 const double t0 = (ft1 * b - ft0 * c) / det;
494 const double t1 = (ft0 * b - ft1 * a) / det;
495
496 if (tvw_ptr)
497 *tvw_ptr = t0;
498 if (tlk_ptr)
499 *tlk_ptr = t1;
500
501 return SOLUTION_EXIST;
502 }
503}
504
506 const double *v_ptr, const int nb, Range edges, double *min_dist_ptr,
507 double *o_ptr, EntityHandle *o_segments) const {
508 MoFEM::Interface &m_field = cOre;
509 moab::Interface &moab(m_field.get_moab());
511
513
514 auto get_point = [i](auto &t_w, auto &t_delta, auto t) {
516 t = std::max(0., std::min(1., t));
517 t_p(i) = t_w(i) + t * t_delta(i);
518 return t_p;
519 };
520
521 auto get_distance = [i](auto &t_p, auto &t_n) {
522 FTensor::Tensor1<double, 3> t_dist_vector;
523 t_dist_vector(i) = t_p(i) - t_n(i);
524 return sqrt(t_dist_vector(i) * t_dist_vector(i));
525 };
526
527 for (auto e : edges) {
528
529 int num_nodes;
530 const EntityHandle *conn_fixed;
531 CHKERR moab.get_connectivity(e, conn_fixed, num_nodes, true);
532 VectorDouble6 coords_fixed(6);
533 CHKERR moab.get_coords(conn_fixed, num_nodes, &coords_fixed[0]);
535 &coords_fixed[0], &coords_fixed[1], &coords_fixed[2]);
537 &coords_fixed[3], &coords_fixed[4], &coords_fixed[5]);
538
539 FTensor::Tensor1<double, 3> t_edge_delta;
540 t_edge_delta(i) = t_f1(i) - t_f0(i);
541
543 v_ptr, v_ptr + 1, v_ptr + 2);
545 o_ptr, o_ptr + 1, o_ptr + 2);
546 FTensor::Tensor0<FTensor::PackPtr<double *, 1>> t_min_dist(min_dist_ptr);
547
548 EntityHandle *colsest_segment_it = nullptr;
549 if (o_segments)
550 colsest_segment_it = o_segments;
551
552 for (int n = 0; n != nb; ++n) {
553
554 double t;
555 if (Tools::minDistancePointFromOnSegment(&t_f0(0), &t_f1(0), &t_n(0),
557 auto t_p = get_point(t_f0, t_edge_delta, t);
558 auto dist_n = get_distance(t_p, t_n);
559 if (dist_n < t_min_dist || t_min_dist < 0) {
560 t_min_dist = dist_n;
561 if (o_ptr)
562 t_min_coords(i) = t_p(i);
563 if (o_segments)
564 *colsest_segment_it = e;
565 }
566 }
567
568 ++t_n;
569 ++t_min_dist;
570 if (o_ptr)
571 ++t_min_coords;
572 if (o_segments)
573 ++colsest_segment_it;
574 }
575 }
576
578}
579
581 MatrixDouble &gauss_pts, const int rule_ksi, const int rule_eta) {
583
584 auto check_rule_edge = [](int rule) {
586 if (rule < 0) {
587 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
588 "Wrong integration rule: %d", rule);
589 }
590 if (rule > QUAD_1D_TABLE_SIZE) {
591 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
592 "rule > quadrature order %d < %d", rule, QUAD_1D_TABLE_SIZE);
593 }
594 if (QUAD_1D_TABLE[rule]->dim != 1) {
595 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong dimension");
596 }
597 if (QUAD_1D_TABLE[rule]->order < rule) {
598 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
599 "wrong order %d != %d", QUAD_1D_TABLE[rule]->order, rule);
600 }
602 };
603
604 CHKERR check_rule_edge(rule_ksi);
605 CHKERR check_rule_edge(rule_eta);
606
607 const int nb_gauss_pts_ksi = QUAD_1D_TABLE[rule_ksi]->npoints;
608 const int nb_gauss_pts_eta = QUAD_1D_TABLE[rule_eta]->npoints;
609 gauss_pts.resize(3, nb_gauss_pts_ksi * nb_gauss_pts_eta, false);
610
611 int gg = 0;
612 for (size_t i = 0; i != nb_gauss_pts_ksi; ++i) {
613 const double wi = QUAD_1D_TABLE[rule_ksi]->weights[i];
614 const double ksi = (QUAD_1D_TABLE[rule_ksi]->points[2 * i + 1]);
615 for (size_t j = 0; j != nb_gauss_pts_eta; ++j, ++gg) {
616 const double wk = wi * QUAD_1D_TABLE[rule_eta]->weights[j];
617 const double eta = QUAD_1D_TABLE[rule_eta]->points[2 * j + 1];
618 gauss_pts(0, gg) = ksi;
619 gauss_pts(1, gg) = eta;
620 gauss_pts(2, gg) = wk;
621 }
622 }
623 if (gg != gauss_pts.size2())
624 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong size of matrix");
625
627}
628
630 MatrixDouble &gauss_pts, const int rule_ksi, const int rule_eta,
631 const int rule_zeta) {
633
634 auto check_rule_edge = [](int rule) {
636 if (rule < 0) {
637 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
638 "Wrong integration rule: %d", rule);
639 }
640 if (rule > QUAD_1D_TABLE_SIZE) {
641 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
642 "rule > quadrature order %d < %d", rule, QUAD_1D_TABLE_SIZE);
643 }
644 if (QUAD_1D_TABLE[rule]->dim != 1) {
645 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong dimension");
646 }
647 if (QUAD_1D_TABLE[rule]->order < rule) {
648 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
649 "wrong order %d != %d", QUAD_1D_TABLE[rule]->order, rule);
650 }
652 };
653
654 CHKERR check_rule_edge(rule_ksi);
655 CHKERR check_rule_edge(rule_eta);
656 CHKERR check_rule_edge(rule_zeta);
657
658 const int nb_gauss_pts_ksi = QUAD_1D_TABLE[rule_ksi]->npoints;
659 const int nb_gauss_pts_eta = QUAD_1D_TABLE[rule_eta]->npoints;
660 const int nb_gauss_pts_zeta = QUAD_1D_TABLE[rule_zeta]->npoints;
661 gauss_pts.resize(4, nb_gauss_pts_ksi * nb_gauss_pts_eta * nb_gauss_pts_zeta,
662 false);
663
664 int gg = 0;
665 for (size_t i = 0; i != nb_gauss_pts_ksi; ++i) {
666 const double wi = QUAD_1D_TABLE[rule_ksi]->weights[i];
667 const double ksi = QUAD_1D_TABLE[rule_ksi]->points[2 * i + 1];
668 for (size_t j = 0; j != nb_gauss_pts_eta; ++j) {
669 const double wj = wi * QUAD_1D_TABLE[rule_eta]->weights[j];
670 const double eta = QUAD_1D_TABLE[rule_eta]->points[2 * j + 1];
671 for (size_t k = 0; k != nb_gauss_pts_zeta; ++k, ++gg) {
672 const double wk = wj * QUAD_1D_TABLE[rule_zeta]->weights[k];
673 const double zeta = QUAD_1D_TABLE[rule_zeta]->points[2 * k + 1];
674 gauss_pts(0, gg) = ksi;
675 gauss_pts(1, gg) = eta;
676 gauss_pts(2, gg) = zeta;
677 gauss_pts(3, gg) = wk;
678 }
679 }
680 }
681
682 if (gg != gauss_pts.size2())
683 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong size of matrix");
684
686}
687
688} // namespace MoFEM
constexpr double a
RowColData
RowColData.
Definition: definitions.h:123
@ COL
Definition: definitions.h:123
@ ROW
Definition: definitions.h:123
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ L2
field with C-1 continuity
Definition: definitions.h:88
#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
Definition: definitions.h:31
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
const int dim
PetscErrorCode ShapeDiffMBTET(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:319
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
PetscErrorCode ShapeFaceNormalMBTRI(double *diffN, const double *coords, double *normal)
Definition: fem_tools.c:229
PetscErrorCode ShapeDiffMBTRI(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:194
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
FTensor::Index< 'n', SPACE_DIM > n
static const double edge_coords[6][6]
constexpr double eta
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
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
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
double tol
FTensor::Tensor1< double, SPACE_DIM > normal(FTensor::Tensor1< T1, 3 > &t_coords, FTensor::Tensor1< T2, SPACE_DIM > &t_disp)
Definition: ContactOps.hpp:65
const FTensor::Tensor2< T, Dim, Dim > Vec
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
auto f
Definition: HenckyOps.hpp:5
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
VectorBoundedArray< double, 6 > VectorDouble6
Definition: Types.hpp:95
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
auto type_from_handle(const EntityHandle h)
get type from entity handle
Definition: Templates.hpp:1479
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1199
constexpr double t
plate stiffness
Definition: plate.cpp:59
static QUAD *const QUAD_1D_TABLE[]
Definition: quad.h:164
#define QUAD_1D_TABLE_SIZE
Definition: quad.h:163
const int N
Definition: speed_test.cpp:3
virtual moab::Interface & get_moab()=0
Deprecated interface functions.
keeps basic data about problem
DofIdx getNbLocalDofsRow() const
DofIdx getNbLocalDofsCol() const
auto & getNumeredColDofsPtr() const
get access to numeredColDofsPtr storing DOFs on cols
auto & getNumeredRowDofsPtr() const
get access to numeredRowDofsPtr storing DOFs on rows
Auxiliary tools.
Definition: Tools.hpp:19
static constexpr double shapeFunMBEDGE1At00
Definition: Tools.hpp:56
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:439
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
static MoFEMErrorCode outerProductOfEdgeIntegrationPtsForQuad(MatrixDouble &pts, const int edge0, const int edge1)
Definition: Tools.cpp:580
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:416
static constexpr double shapeFunMBTRI1At00
Definition: Tools.hpp:113
static constexpr double shapeFunMBTRI0At00
Definition: Tools.hpp:112
static constexpr std::array< double, 8 > diffShapeFunMBQUADAtCenter
Definition: Tools.hpp:212
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:505
static constexpr std::array< double, 6 > diffShapeFunMBTRI
Definition: Tools.hpp:104
SEGMENT_MIN_DISTANCE
Definition: Tools.hpp:507
@ SOLUTION_EXIST
Definition: Tools.hpp:508
@ SEGMENT_TWO_AND_TWO_ARE_POINT
Definition: Tools.hpp:511
@ SEGMENT_ONE_IS_POINT
Definition: Tools.hpp:509
@ SEGMENT_TWO_IS_POINT
Definition: Tools.hpp:510
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
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
static constexpr std::array< double, 2 > shapeFunMBEDGEAt00
Array of shape function at zero local point on reference element.
Definition: Tools.hpp:62
static constexpr double shapeFunMBEDGE0At00
Definition: Tools.hpp:55
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
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
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
static MoFEMErrorCode outerProductOfEdgeIntegrationPtsForHex(MatrixDouble &pts, const int edge0, const int edge1, const int edge2)
Definition: Tools.cpp:629
static constexpr double shapeFunMBTRI2At00
Definition: Tools.hpp:114
static constexpr std::array< double, 24 > diffShapeFunMBHEXAtCenter
Definition: Tools.hpp:218
static constexpr std::array< double, 4 > shapeFunMBTETAt000
Array of shape function at zero local point on reference element.
Definition: Tools.hpp:330
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
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
static constexpr std::array< double, 2 > diffShapeFunMBEDGE
Definition: Tools.hpp:68
static double tetVolume(const double *coords)
Calculate volume of tetrahedron.
Definition: Tools.cpp:32
static double volumeLengthQuality(const double *coords)
Calculate tetrahedron volume length quality.
Definition: Tools.cpp:17
double getTriArea(const EntityHandle tri) const
Get triangle area.
Definition: Tools.cpp:377
static constexpr std::array< double, 3 > shapeFunMBTRIAt00
Array of shape function at zero local point on reference element.
Definition: Tools.hpp:120
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
Definition: Tools.cpp:11
static constexpr std::array< double, 12 > diffShapeFunMBTET
Definition: Tools.hpp:271
static double getEdgeLength(const double *edge_coords)
Get edge length.
Definition: Tools.cpp:385
MoFEM::Core & cOre
Definition: Tools.hpp:24
static MoFEMErrorCode getTriNormal(const double *coords, double *normal)
Get the Tri Normal objectGet triangle normal.
Definition: Tools.cpp:352
base class for all interface classes
double * points
Definition: quad.h:30
double * weights
Definition: quad.h:31
int npoints
Definition: quad.h:29