v0.9.0
Tools.cpp
Go to the documentation of this file.
1 /** \file Tools.cpp
2  * \brief Auxiliary tools
3  *
4  * MoFEM is free software: you can redistribute it and/or modify it under
5  * the terms of the GNU Lesser General Public License as published by the
6  * Free Software Foundation, either version 3 of the License, or (at your
7  * option) any later version.
8  *
9  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
11  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
12  * License for more details.
13  *
14  * You should have received a copy of the GNU Lesser General Public
15  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>
16  */
17 
18 #include <phg-quadrule/quad.h>
19 
20 namespace MoFEM {
21 
23  UnknownInterface **iface) const {
25  *iface = NULL;
26  if (uuid == IDD_MOFEMNodeMerger) {
27  *iface = const_cast<Tools *>(this);
29  }
30  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "unknown interface");
32 }
33 
34 double Tools::volumeLengthQuality(const double *coords) {
35  double lrms = 0;
36  for (int dd = 0; dd != 3; dd++) {
37  lrms += pow(coords[0 * 3 + dd] - coords[1 * 3 + dd], 2) +
38  pow(coords[0 * 3 + dd] - coords[2 * 3 + dd], 2) +
39  pow(coords[0 * 3 + dd] - coords[3 * 3 + dd], 2) +
40  pow(coords[1 * 3 + dd] - coords[2 * 3 + dd], 2) +
41  pow(coords[1 * 3 + dd] - coords[3 * 3 + dd], 2) +
42  pow(coords[2 * 3 + dd] - coords[3 * 3 + dd], 2);
43  }
44  lrms = sqrt((1. / 6.) * lrms);
45  double volume = tetVolume(coords);
46  return 6. * sqrt(2.) * volume / pow(lrms, 3);
47 }
48 
49 double Tools::tetVolume(const double *coords) {
50  double diff_n[12];
51  ShapeDiffMBTET(diff_n);
52  FTensor::Tensor1<double *, 3> t_diff_n(&diff_n[0], &diff_n[1], &diff_n[2], 3);
53  FTensor::Tensor1<const double *, 3> t_coords(&coords[0], &coords[1],
54  &coords[2], 3);
58  jac(i, j) = 0;
59  for (int nn = 0; nn != 4; nn++) {
60  jac(i, j) += t_coords(i) * t_diff_n(j);
61  ++t_coords;
62  ++t_diff_n;
63  }
64  return dEterminant(jac) / 6.;
65 }
66 
68 Tools::minTetsQuality(const Range &tets, double &min_quality, Tag th,
69  boost::function<double(double, double)> f) {
70  MoFEM::Interface &m_field = cOre;
71  moab::Interface &moab(m_field.get_moab());
73  const EntityHandle *conn;
74  int num_nodes;
75  double coords[12];
76  for (auto tet : tets) {
77  CHKERR m_field.get_moab().get_connectivity(tet, conn, num_nodes, true);
78  if (th) {
79  CHKERR moab.tag_get_data(th, conn, num_nodes, coords);
80  } else {
81  CHKERR moab.get_coords(conn, num_nodes, coords);
82  }
83  double q = Tools::volumeLengthQuality(coords);
84  if (!std::isnormal(q))
85  q = -2;
86  min_quality = f(q, min_quality);
87  }
89 }
90 
91 constexpr std::array<double, 2> Tools::diffShapeFunMBEDGE;
92 constexpr std::array<double, 6> Tools::diffShapeFunMBTRI;
93 constexpr std::array<double, 12> Tools::diffShapeFunMBTET;
94 constexpr std::array<double, 4> Tools::shapeFunMBTETAt000;
95 constexpr std::array<double, 8> Tools::diffShapeFunMBQUADAtCenter;
96 
98  const double *elem_coords, const double *global_coords, const int nb_nodes,
99  double *local_coords) {
102 
104  &elem_coords[0], &elem_coords[3], &elem_coords[6], &elem_coords[9]};
105 
108  shapeFunMBTETAt000[3]};
109  FTensor::Tensor1<double, 3> t_coords_at_0;
110 
111  // Build matrix and get coordinates of zero point
112  // ii - global coordinates
113  // jj - local direvatives
114  MatrixDouble3by3 a(3, 3);
115  for (auto ii : {0, 1, 2}) {
118  &diffShapeFunMBTET[9]);
119  for (auto jj : {0, 1, 2}) {
120  a(jj, ii) = t_diff(i) * t_elem_coords(i);
121  ++t_diff;
122  }
123  t_coords_at_0(ii) = t_n(i) * t_elem_coords(i);
124  ++t_elem_coords;
125  }
126 
128  &global_coords[0], &global_coords[1], &global_coords[2]};
130  &local_coords[0], &local_coords[1], &local_coords[2]};
131 
132  // Calculate right hand side
134  for (int ii = 0; ii != nb_nodes; ++ii) {
135  t_local_coords(j) = t_global_coords(j) - t_coords_at_0(j);
136  ++t_local_coords;
137  ++t_global_coords;
138  }
139 
140  // Solve problem
141  int IPIV[3];
142  int info = lapack_dgesv(3, nb_nodes, &a(0, 0), 3, IPIV, local_coords, 3);
143  if (info != 0)
144  SETERRQ1(PETSC_COMM_SELF, 1, "info == %d", info);
145 
147 }
148 
149 MoFEMErrorCode Tools::getTetsWithQuality(Range &out_tets, const Range &tets,
150  Tag th,
151  boost::function<bool(double)> f) {
152  MoFEM::Interface &m_field = cOre;
153  moab::Interface &moab(m_field.get_moab());
155  Range to_write;
156  const EntityHandle *conn;
157  int num_nodes;
158  double coords[12];
159  for (auto tet : tets) {
160  CHKERR m_field.get_moab().get_connectivity(tet, conn, num_nodes, true);
161  if (th) {
162  CHKERR moab.tag_get_data(th, conn, num_nodes, coords);
163  } else {
164  CHKERR moab.get_coords(conn, num_nodes, coords);
165  }
166  double q = Tools::volumeLengthQuality(coords);
167  if (f(q)) {
168  out_tets.insert(tet);
169  }
170  }
172 }
173 
175  const char *file_type,
176  const char *options,
177  const Range &tets, Tag th,
178  boost::function<bool(double)> f) {
179  MoFEM::Interface &m_field = cOre;
180  moab::Interface &moab(m_field.get_moab());
182  Range out_tets;
183  CHKERR getTetsWithQuality(out_tets, tets, th, f);
184  EntityHandle meshset;
185  CHKERR moab.create_meshset(MESHSET_SET, meshset);
186  CHKERR moab.add_entities(meshset, out_tets);
187  CHKERR moab.write_file(file_name, file_type, options, &meshset, 1);
188  CHKERR moab.delete_entities(&meshset, 1);
190 }
191 
192 MoFEMErrorCode Tools::checkIfPointIsInTet(const double tet_coords[],
193  const double global_coord[],
194  const double tol, bool &result) {
195  double loc_coord[] = {0, 0, 0};
196  double N[4], diffN[12];
198  CHKERR ShapeDiffMBTET(diffN);
199  CHKERR ShapeMBTET(N, &loc_coord[0], &loc_coord[1], &loc_coord[2], 1);
200  CHKERR ShapeMBTET_inverse(N, diffN, tet_coords, global_coord, loc_coord);
201  CHKERR ShapeMBTET(N, &loc_coord[0], &loc_coord[1], &loc_coord[2], 1);
202  result = true;
203  for (int n = 0; n != 4; ++n) {
204  if (N[n] < -tol || (N[n] - 1) > tol) {
205  result = false;
206  break;
207  }
208  }
210 }
211 
213  const RowColData row_or_col,
214  Vec v) {
216  int loc_size;
217  CHKERR VecGetLocalSize(v, &loc_size);
218  int prb_loc_size = 0;
219  boost::shared_ptr<NumeredDofEntity_multiIndex> prb_dofs;
220  switch (row_or_col) {
221  case ROW:
222  prb_loc_size = prb_ptr->getNbLocalDofsRow();
223  prb_dofs = prb_ptr->getNumeredDofsRows();
224  break;
225  case COL:
226  prb_loc_size = prb_ptr->getNbLocalDofsCol();
227  prb_dofs = prb_ptr->getNumeredDofsCols();
228  break;
229  break;
230  default:
231  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
232  "Wrong argument, row_or_col should be row or column");
233  }
234  if (loc_size != prb_loc_size) {
235  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
236  "Inconsistent size of vector and problem %d != %d", loc_size,
237  prb_loc_size);
238  }
239  const double *a;
240  CHKERR VecGetArrayRead(v, &a);
241  MPI_Comm comm = PetscObjectComm((PetscObject)v);
242  for (int ii = 0; ii != loc_size; ++ii) {
243  if (!boost::math::isfinite(a[ii])) {
244  NumeredDofEntityByLocalIdx::iterator dit =
245  prb_dofs->get<PetscLocalIdx_mi_tag>().find(ii);
246  std::ostringstream ss;
247  ss << "Not a number " << a[ii] << " on dof: " << endl
248  << **dit << endl
249  << endl;
250  PetscSynchronizedPrintf(comm, "%s", ss.str().c_str());
251  }
252  }
253  CHKERR VecRestoreArrayRead(v, &a);
254  PetscSynchronizedFlush(comm, PETSC_STDOUT);
256 }
257 
258 MoFEMErrorCode Tools::getTriNormal(const double *coords, double *normal) {
260  double diffN[6];
261  CHKERR ShapeDiffMBTRI(diffN);
262  CHKERR ShapeFaceNormalMBTRI(diffN, coords, normal);
264 }
265 
267  double *normal) const {
268  MoFEM::Interface &m_field = cOre;
269  moab::Interface &moab(m_field.get_moab());
271  if (moab.type_from_handle(tri) != MBTRI) {
272  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "Works only for triangle");
273  }
274  const EntityHandle *conn;
275  int num_nodes;
276  double coords[9];
277  CHKERR moab.get_connectivity(tri, conn, num_nodes, true);
278  CHKERR moab.get_coords(conn, num_nodes, coords);
279  CHKERR getTriNormal(coords, normal);
281 }
282 
283 double Tools::getTriArea(const EntityHandle tri) const {
285  ierr = getTriNormal(tri, &t_normal(0));
286  CHKERRABORT(PETSC_COMM_SELF, ierr);
288  return sqrt(t_normal(i) * t_normal(i)) * 0.5;
289 }
290 
291 double Tools::getEdgeLength(const double *edge_coords) {
293  edge_coords[2]);
295  edge_coords[5]);
297  t_coords_n0(i) -= t_coords_n1(i);
298  return sqrt(t_coords_n0(i) * t_coords_n0(i));
299 }
300 
301 double Tools::getEdgeLength(const EntityHandle edge) {
302  MoFEM::Interface &m_field = cOre;
303  moab::Interface &moab(m_field.get_moab());
304  auto get_edge_coords = [edge, &moab](double *const coords) {
306  if (moab.type_from_handle(edge) != MBEDGE) {
307  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA, "Works only for edge");
308  }
309  const EntityHandle *conn;
310  int num_nodes;
311  CHKERR moab.get_connectivity(edge, conn, num_nodes, true);
312  CHKERR moab.get_coords(conn, 2, coords);
314  };
315  double coords[6];
316  ierr = get_edge_coords(coords);
317  CHKERRABORT(PETSC_COMM_SELF, ierr);
318  return getEdgeLength(coords);
319 }
320 
322 Tools::minDistancePointFromOnSegment(const double *w_ptr, const double *v_ptr,
323  const double *p_ptr, double *const t_ptr) {
325  FTensor::Tensor1<const double *, 3> t_w(&w_ptr[0], &w_ptr[1], &w_ptr[2]);
326  FTensor::Tensor1<const double *, 3> t_v(&v_ptr[0], &v_ptr[1], &v_ptr[2]);
327  FTensor::Tensor1<const double *, 3> t_p(&p_ptr[0], &p_ptr[1], &p_ptr[2]);
329  t_vw(i) = t_v(i) - t_w(i);
330  const double dot_vw = t_vw(i) * t_vw(i);
331  if (dot_vw == 0) {
332  if (t_ptr)
333  *t_ptr = 0;
334  return SEGMENT_ONE_IS_POINT;
335  }
337  t_pw(i) = t_p(i) - t_w(i);
338  const double t = t_pw(i) * t_vw(i) / dot_vw;
339  if (t_ptr)
340  *t_ptr = t;
341  return SOLUTION_EXIST;
342 }
343 
345 Tools::minDistanceFromSegments(const double *w_ptr, const double *v_ptr,
346  const double *k_ptr, const double *l_ptr,
347  double *const tvw_ptr, double *const tlk_ptr) {
348 
350  FTensor::Tensor1<const double *, 3> t_w(&w_ptr[0], &w_ptr[1], &w_ptr[2]);
351  FTensor::Tensor1<const double *, 3> t_v(&v_ptr[0], &v_ptr[1], &v_ptr[2]);
352  FTensor::Tensor1<const double *, 3> t_k(&k_ptr[0], &k_ptr[1], &k_ptr[2]);
353  FTensor::Tensor1<const double *, 3> t_l(&l_ptr[0], &l_ptr[1], &l_ptr[2]);
354 
355  // First segnent is a point
357  t_vw(i) = t_v(i) - t_w(i);
358  double dot_vw = t_vw(i) * t_vw(i);
359  if (dot_vw == 0) {
360  if (tvw_ptr)
361  *tvw_ptr = 0;
362  if (minDistancePointFromOnSegment(k_ptr, l_ptr, w_ptr, tlk_ptr) ==
365  else
366  return SEGMENT_ONE_IS_POINT;
367  }
368 
369  // Second segment is a point
371  t_lk(i) = t_l(i) - t_k(i);
372  double dot_lk = t_lk(i) * t_lk(i);
373  if (dot_lk == 0) {
374  if (tlk_ptr)
375  *tlk_ptr = 0;
376  if (minDistancePointFromOnSegment(w_ptr, v_ptr, k_ptr, tvw_ptr) ==
379  else
380  return SEGMENT_TWO_IS_POINT;
381  }
382 
383  const double a = t_vw(i) * t_vw(i);
384  const double b = -t_vw(i) * t_lk(i);
385  const double c = t_lk(i) * t_lk(i);
386 
387  const double det = a * c - b * b;
388  if (det == 0) {
389 
390  return NO_SOLUTION;
391 
392  } else {
393 
395  t_wk(i) = t_w(i) - t_k(i);
396 
397  const double ft0 = t_vw(i) * t_wk(i);
398  const double ft1 = -t_lk(i) * t_wk(i);
399  const double t0 = (ft1 * b - ft0 * c) / det;
400  const double t1 = (ft0 * b - ft1 * a) / det;
401 
402  if (tvw_ptr)
403  *tvw_ptr = t0;
404  if (tlk_ptr)
405  *tlk_ptr = t1;
406 
407  return SOLUTION_EXIST;
408  }
409 }
410 
412  const double *v_ptr, const int nb, Range edges, double *min_dist_ptr,
413  double *o_ptr, EntityHandle *o_segments) const {
414  MoFEM::Interface &m_field = cOre;
415  moab::Interface &moab(m_field.get_moab());
417 
419 
420  auto get_point = [i](auto &t_w, auto &t_delta, auto t) {
422  t = std::max(0., std::min(1., t));
423  t_p(i) = t_w(i) + t * t_delta(i);
424  return t_p;
425  };
426 
427  auto get_distance = [i](auto &t_p, auto &t_n) {
428  FTensor::Tensor1<double, 3> t_dist_vector;
429  t_dist_vector(i) = t_p(i) - t_n(i);
430  return sqrt(t_dist_vector(i) * t_dist_vector(i));
431  };
432 
433  for (auto e : edges) {
434 
435  int num_nodes;
436  const EntityHandle *conn_fixed;
437  CHKERR moab.get_connectivity(e, conn_fixed, num_nodes, true);
438  VectorDouble6 coords_fixed(6);
439  CHKERR moab.get_coords(conn_fixed, num_nodes, &coords_fixed[0]);
441  &coords_fixed[0], &coords_fixed[1], &coords_fixed[2]);
443  &coords_fixed[3], &coords_fixed[4], &coords_fixed[5]);
444 
445  FTensor::Tensor1<double, 3> t_edge_delta;
446  t_edge_delta(i) = t_f1(i) - t_f0(i);
447 
449  v_ptr, v_ptr + 1, v_ptr + 2);
451  o_ptr, o_ptr + 1, o_ptr + 2);
452  FTensor::Tensor0<FTensor::PackPtr<double *, 1>> t_min_dist(min_dist_ptr);
453 
454  EntityHandle *colsest_segment_it = nullptr;
455  if (o_segments)
456  colsest_segment_it = o_segments;
457 
458  for (int n = 0; n != nb; ++n) {
459 
460  double t;
461  if (Tools::minDistancePointFromOnSegment(&t_f0(0), &t_f1(0), &t_n(0),
462  &t) == Tools::SOLUTION_EXIST) {
463  auto t_p = get_point(t_f0, t_edge_delta, t);
464  auto dist_n = get_distance(t_p, t_n);
465  if (dist_n < t_min_dist || t_min_dist < 0) {
466  t_min_dist = dist_n;
467  if (o_ptr)
468  t_min_coords(i) = t_p(i);
469  if (o_segments)
470  *colsest_segment_it = e;
471  }
472  }
473 
474  ++t_n;
475  ++t_min_dist;
476  if (o_ptr)
477  ++t_min_coords;
478  if (o_segments)
479  ++colsest_segment_it;
480  }
481  }
482 
484 }
485 
487  MatrixDouble &gauss_pts, const int rule_ksi, const int rule_eta) {
489 
490  auto check_rule_edge = [](int rule) {
492  if (rule < 0) {
493  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
494  "Wrong integration rule: %d", rule);
495  }
496  if (rule > QUAD_1D_TABLE_SIZE) {
497  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
498  "rule > quadrature order %d < %d", rule, QUAD_1D_TABLE_SIZE);
499  }
500  if (QUAD_1D_TABLE[rule]->dim != 1) {
501  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong dimension");
502  }
503  if (QUAD_1D_TABLE[rule]->order < rule) {
504  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
505  "wrong order %d != %d", QUAD_1D_TABLE[rule]->order, rule);
506  }
508  };
509 
510  CHKERR check_rule_edge(rule_ksi);
511  CHKERR check_rule_edge(rule_eta);
512 
513  int nb_gauss_pts_ksi = QUAD_1D_TABLE[rule_ksi]->npoints;
514  int nb_gauss_pts_eta = QUAD_1D_TABLE[rule_eta]->npoints;
515  gauss_pts.resize(3, nb_gauss_pts_ksi * nb_gauss_pts_eta, false);
516  gauss_pts.clear();
517 
518  VectorDouble ones;
519  ones.resize(max(nb_gauss_pts_ksi, nb_gauss_pts_eta), false);
520  fill(ones.begin(), ones.end(), 1.0);
521 
522  cblas_dger(CblasRowMajor, nb_gauss_pts_eta, nb_gauss_pts_ksi, 1, &ones(0), 1,
523  &QUAD_1D_TABLE[rule_ksi]->points[1], 2, &gauss_pts(0, 0),
524  nb_gauss_pts_ksi);
525 
526  cblas_dger(CblasRowMajor, nb_gauss_pts_eta, nb_gauss_pts_ksi, 1,
527  &QUAD_1D_TABLE[rule_eta]->points[1], 2, &ones(0), 1,
528  &gauss_pts(1, 0), nb_gauss_pts_ksi);
529 
530  cblas_dger(CblasRowMajor, nb_gauss_pts_eta, nb_gauss_pts_ksi, 1,
531  QUAD_1D_TABLE[rule_eta]->weights, 1,
532  QUAD_1D_TABLE[rule_ksi]->weights, 1, &gauss_pts(2, 0),
533  nb_gauss_pts_ksi);
534 
536 }
537 
538 } // namespace MoFEM
MoFEM interface unique ID.
static double volumeLengthQuality(const double *coords)
Calculate tetrahedron volume length quality.
Definition: Tools.cpp:34
static QUAD *const QUAD_1D_TABLE[]
Definition: quad.h:164
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:411
virtual moab::Interface & get_moab()=0
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:318
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:149
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
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:345
int npoints
Definition: quad.h:29
PetscErrorCode ShapeFaceNormalMBTRI(double *diffN, const double *coords, double *normal)
Definition: fem_tools.c:241
static constexpr std::array< double, 2 > diffShapeFunMBEDGE
Definition: Tools.hpp:71
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
base class for all interface classes
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
static constexpr std::array< double, 12 > diffShapeFunMBTET
Definition: Tools.hpp:143
double getTriArea(const EntityHandle tri) const
Get triangle area.
Definition: Tools.cpp:283
DofIdx getNbLocalDofsRow() const
static const MOFEMuuid IDD_MOFEMNodeMerger
Definition: NodeMerger.hpp:23
void cblas_dger(const enum CBLAS_ORDER order, const int M, const int N, const double alpha, const double *X, const int incX, const double *Y, const int incY, double *A, const int lda)
Definition: cblas_dger.c:12
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:97
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
keeps basic data about problemThis is low level structure with information about problem,...
RowColData
RowColData.
Definition: definitions.h:186
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:174
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
const boost::shared_ptr< NumeredDofEntity_multiIndex > & getNumeredDofsCols() const
get access to numeredDofsCols storing DOFs on cols
static constexpr std::array< double, 4 > shapeFunMBTETAt000
Array of shape function at zero local point on reference element.
Definition: Tools.hpp:202
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:192
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:347
static double tetVolume(const double *coords)
Calculate volume of tetrahedron.
Definition: Tools.cpp:49
VectorBoundedArray< double, 6 > VectorDouble6
Definition: Types.hpp:88
PetscErrorCode ShapeDiffMBTET(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:331
PetscErrorCode ShapeDiffMBTRI(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:206
double tol
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:68
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
static double dEterminant(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:385
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
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:212
static MoFEMErrorCode getTriNormal(const double *coords, double *normal)
Get the Tri Normal objectGet triangle normal.
Definition: Tools.cpp:258
MoFEMErrorCode query_interface(const MOFEMuuid &uuid, UnknownInterface **iface) const
Definition: Tools.cpp:22
static const double edge_coords[6][6]
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition: Types.hpp:95
static constexpr std::array< double, 8 > diffShapeFunMBQUADAtCenter
Definition: Tools.hpp:112
static constexpr std::array< double, 6 > diffShapeFunMBTRI
Definition: Tools.hpp:87
#define CHKERR
Inline error check.
Definition: definitions.h:596
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:322
static MoFEMErrorCode outerProductOfEdgeIntegrationPtsForQuad(MatrixDouble &pts, const int edge0, const int edge1)
Definition: Tools.cpp:486
static double getEdgeLength(const double *edge_coords)
Get edge length.
Definition: Tools.cpp:291
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:72
DofIdx getNbLocalDofsCol() const
const int N
Definition: speed_test.cpp:3
MoFEM::Core & cOre
Definition: Tools.hpp:35
const boost::shared_ptr< NumeredDofEntity_multiIndex > & getNumeredDofsRows() const
get access to numeredDofsRows storing DOFs on rows
static __CLPK_integer lapack_dgesv(__CLPK_integer n, __CLPK_integer nrhs, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *b, __CLPK_integer ldb)
Definition: lapack_wrap.h:188
SEGMENT_MIN_DISTANCE
Definition: Tools.hpp:332
#define QUAD_1D_TABLE_SIZE
Definition: quad.h:163