13 *iface =
const_cast<Tools *
>(
this);
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);
27 lrms = sqrt((1. / 6.) * lrms);
29 return 6. * sqrt(2.) * volume / pow(lrms, 3);
42 for (
int nn = 0; nn != 4; nn++) {
43 jac(
i,
j) += t_coords(
i) * t_diff_n(
j);
52 boost::function<
double(
double,
double)> f) {
54 moab::Interface &moab(m_field.
get_moab());
59 for (
auto tet : tets) {
60 CHKERR m_field.
get_moab().get_connectivity(tet, conn, num_nodes,
true);
62 CHKERR moab.tag_get_data(
th, conn, num_nodes, coords);
64 CHKERR moab.get_coords(conn, num_nodes, coords);
67 if (!std::isnormal(q))
69 min_quality = f(q, min_quality);
91 const double *elem_coords,
const double *global_coords,
const int nb_nodes,
92 double *local_coords) {
97 &elem_coords[0], &elem_coords[3], &elem_coords[6], &elem_coords[9]};
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);
116 t_coords_at_0(ii) = t_n(
i) * t_elem_coords(
i);
121 &global_coords[0], &global_coords[1], &global_coords[2]};
123 &local_coords[0], &local_coords[1], &local_coords[2]};
127 for (
int ii = 0; ii != nb_nodes; ++ii) {
128 t_local_coords(
j) = t_global_coords(
j) - t_coords_at_0(
j);
135 int info =
lapack_dgesv(3, nb_nodes, &
a(0, 0), 3, IPIV, local_coords, 3);
137 SETERRQ1(PETSC_COMM_SELF, 1,
"info == %d", info);
143 const double *elem_coords,
const double *global_coords,
const int nb_nodes,
144 double *local_coords) {
153 &elem_coords[0], &elem_coords[3], &elem_coords[6]};
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);
170 t_coords_at_0(ii) = t_n(i3) * t_elem_coords(i3);
175 &global_coords[0], &global_coords[1], &global_coords[2]};
177 &local_coords[0], &local_coords[1]};
180 t_b(K2,
L2) = t_a(K2, j3) ^ t_a(
L2, j3);
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;
189 for (
int ii = 0; ii != nb_nodes; ++ii) {
192 (t_a(K2, j3) * (t_global_coords(j3) - t_coords_at_0(j3)));
201 const double *elem_coords,
const double *global_coords,
const int nb_nodes,
202 double *local_coords) {
208 &elem_coords[0], &elem_coords[3], &elem_coords[6]};
216 for (
auto ii : {0, 1, 2}) {
225 &global_coords[0], &global_coords[1], &global_coords[2]};
229 const double b = t_a(i3) * t_a(i3);
230 const double inv_b = 1 / b;
233 for (
int ii = 0; ii != nb_nodes; ++ii) {
235 inv_b * (t_a(i3) * (t_global_coords(i3) - t_coords_at_0(i3)));
245 boost::function<
bool(
double)> f) {
247 moab::Interface &moab(m_field.
get_moab());
253 for (
auto tet : tets) {
254 CHKERR m_field.
get_moab().get_connectivity(tet, conn, num_nodes,
true);
256 CHKERR moab.tag_get_data(
th, conn, num_nodes, coords);
258 CHKERR moab.get_coords(conn, num_nodes, coords);
262 out_tets.insert(tet);
269 const char *file_type,
272 boost::function<
bool(
double)> f) {
274 moab::Interface &moab(m_field.
get_moab());
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);
287 const double global_coord[],
288 const double tol,
bool &result) {
289 double loc_coord[] = {0, 0, 0};
290 double N[4], diffN[12];
297 for (
int n = 0;
n != 4; ++
n) {
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) {
326 "Wrong argument, row_or_col should be row or column");
328 if (loc_size != prb_loc_size) {
330 "Inconsistent size of vector and problem %d != %d", loc_size,
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 =
340 std::ostringstream ss;
341 ss <<
"Not a number " <<
a[ii] <<
" on dof: " << endl
344 PetscSynchronizedPrintf(comm,
"%s", ss.str().c_str());
348 PetscSynchronizedFlush(comm, PETSC_STDOUT);
361 double *normal)
const {
363 moab::Interface &moab(m_field.
get_moab());
371 CHKERR moab.get_connectivity(tri, conn, num_nodes,
true);
372 CHKERR moab.get_coords(conn, num_nodes, coords);
381 return sqrt(t_normal(
i) * t_normal(
i)) * 0.5;
390 t_coords_n0(
i) -= t_coords_n1(
i);
391 return sqrt(t_coords_n0(
i) * t_coords_n0(
i));
396 moab::Interface &moab(m_field.
get_moab());
397 auto get_edge_coords = [edge, &moab](
double *
const coords) {
404 CHKERR moab.get_connectivity(edge, conn, num_nodes,
true);
405 CHKERR moab.get_coords(conn, 2, coords);
409 ierr = get_edge_coords(coords);
410 CHKERRABORT(PETSC_COMM_SELF,
ierr);
416 const double *p_ptr,
double *
const t_ptr) {
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()) {
430 t_pw(
i) = t_p(
i) - t_w(
i);
431 const double t = t_pw(
i) * t_vw(
i) / dot_vw;
439 const double *k_ptr,
const double *l_ptr,
440 double *
const tvw_ptr,
double *
const tlk_ptr) {
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()) {
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()) {
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);
480 const double det =
a *
c - b * b;
481 if (std::abs(det) < std::numeric_limits<double>::epsilon()) {
488 t_wk(
i) = t_w(
i) - t_k(
i);
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;
505 const double *v_ptr,
const int nb,
Range edges,
double *min_dist_ptr,
508 moab::Interface &moab(m_field.
get_moab());
513 auto get_point = [
i](
auto &t_w,
auto &t_delta,
auto t) {
515 t = std::max(0., std::min(1.,
t));
516 t_p(
i) = t_w(
i) +
t * t_delta(
i);
520 auto get_distance = [
i](
auto &t_p,
auto &t_n) {
522 t_dist_vector(
i) = t_p(
i) - t_n(
i);
523 return sqrt(t_dist_vector(
i) * t_dist_vector(
i));
526 for (
auto e : edges) {
530 CHKERR moab.get_connectivity(e, conn_fixed, num_nodes,
true);
532 CHKERR moab.get_coords(conn_fixed, num_nodes, &coords_fixed[0]);
534 &coords_fixed[0], &coords_fixed[1], &coords_fixed[2]);
536 &coords_fixed[3], &coords_fixed[4], &coords_fixed[5]);
539 t_edge_delta(
i) = t_f1(
i) - t_f0(
i);
542 v_ptr, v_ptr + 1, v_ptr + 2);
544 o_ptr, o_ptr + 1, o_ptr + 2);
549 colsest_segment_it = o_segments;
551 for (
int n = 0;
n != nb; ++
n) {
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) {
561 t_min_coords(
i) = t_p(
i);
563 *colsest_segment_it = e;
572 ++colsest_segment_it;
580 MatrixDouble &gauss_pts,
const int rule_ksi,
const int rule_eta) {
583 auto check_rule_edge = [](
int rule) {
587 "Wrong integration rule: %d", rule);
603 CHKERR check_rule_edge(rule_ksi);
604 CHKERR check_rule_edge(rule_eta);
608 gauss_pts.resize(3, nb_gauss_pts_ksi * nb_gauss_pts_eta,
false);
611 for (
size_t i = 0;
i != nb_gauss_pts_ksi; ++
i) {
614 for (
size_t j = 0;
j != nb_gauss_pts_eta; ++
j, ++gg) {
617 gauss_pts(0, gg) = ksi;
618 gauss_pts(1, gg) =
eta;
619 gauss_pts(2, gg) = wk;
622 if (gg != gauss_pts.size2())
629 MatrixDouble &gauss_pts,
const int rule_ksi,
const int rule_eta,
630 const int rule_zeta) {
633 auto check_rule_edge = [](
int rule) {
637 "Wrong integration rule: %d", rule);
653 CHKERR check_rule_edge(rule_ksi);
654 CHKERR check_rule_edge(rule_eta);
655 CHKERR check_rule_edge(rule_zeta);
660 gauss_pts.resize(4, nb_gauss_pts_ksi * nb_gauss_pts_eta * nb_gauss_pts_zeta,
664 for (
size_t i = 0;
i != nb_gauss_pts_ksi; ++
i) {
667 for (
size_t j = 0;
j != nb_gauss_pts_eta; ++
j) {
670 for (
size_t k = 0;
k != nb_gauss_pts_zeta; ++
k, ++gg) {
673 gauss_pts(0, gg) = ksi;
674 gauss_pts(1, gg) =
eta;
675 gauss_pts(2, gg) =
zeta;
676 gauss_pts(3, gg) = wk;
681 if (gg != gauss_pts.size2())
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ L2
field with C-1 continuity
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
FTensor::Index< 'n', SPACE_DIM > n
static const double edge_coords[6][6]
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
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)
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
VectorBoundedArray< double, 6 > VectorDouble6
implementation of Data Operators for Forces and Sources
auto type_from_handle(const EntityHandle h)
get type from entity handle
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
constexpr double t
plate stiffness
static QUAD *const QUAD_1D_TABLE[]
#define QUAD_1D_TABLE_SIZE
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
base class for all interface classes