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]};
214 for (
auto ii : {0, 1, 2}) {
223 &global_coords[0], &global_coords[1], &global_coords[2]};
227 const double b = t_a(i3) * t_a(i3);
228 const double inv_b = 1 / b;
231 for (
int ii = 0; ii != nb_nodes; ++ii) {
233 inv_b * (t_a(i3) * (t_global_coords(i3) - t_coords_at_0(i3)));
243 boost::function<
bool(
double)> f) {
245 moab::Interface &moab(m_field.
get_moab());
251 for (
auto tet : tets) {
252 CHKERR m_field.
get_moab().get_connectivity(tet, conn, num_nodes,
true);
254 CHKERR moab.tag_get_data(
th, conn, num_nodes, coords);
256 CHKERR moab.get_coords(conn, num_nodes, coords);
260 out_tets.insert(tet);
267 const char *file_type,
270 boost::function<
bool(
double)> f) {
272 moab::Interface &moab(m_field.
get_moab());
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);
285 const double global_coord[],
286 const double tol,
bool &result) {
287 double loc_coord[] = {0, 0, 0};
288 double N[4], diffN[12];
295 for (
int n = 0;
n != 4; ++
n) {
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) {
324 "Wrong argument, row_or_col should be row or column");
326 if (loc_size != prb_loc_size) {
328 "Inconsistent size of vector and problem %d != %d", loc_size,
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 =
338 std::ostringstream ss;
339 ss <<
"Not a number " <<
a[ii] <<
" on dof: " << endl
342 PetscSynchronizedPrintf(comm,
"%s", ss.str().c_str());
346 PetscSynchronizedFlush(comm, PETSC_STDOUT);
359 double *normal)
const {
361 moab::Interface &moab(m_field.
get_moab());
369 CHKERR moab.get_connectivity(tri, conn, num_nodes,
true);
370 CHKERR moab.get_coords(conn, num_nodes, coords);
379 return sqrt(t_normal(
i) * t_normal(
i)) * 0.5;
388 t_coords_n0(
i) -= t_coords_n1(
i);
389 return sqrt(t_coords_n0(
i) * t_coords_n0(
i));
394 moab::Interface &moab(m_field.
get_moab());
395 auto get_edge_coords = [edge, &moab](
double *
const coords) {
402 CHKERR moab.get_connectivity(edge, conn, num_nodes,
true);
403 CHKERR moab.get_coords(conn, 2, coords);
407 ierr = get_edge_coords(coords);
408 CHKERRABORT(PETSC_COMM_SELF,
ierr);
414 const double *p_ptr,
double *
const t_ptr) {
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()) {
428 t_pw(
i) = t_p(
i) - t_w(
i);
429 const double t = t_pw(
i) * t_vw(
i) / dot_vw;
437 const double *k_ptr,
const double *l_ptr,
438 double *
const tvw_ptr,
double *
const tlk_ptr) {
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()) {
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()) {
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);
478 const double det =
a *
c - b * b;
479 if (std::abs(det) < std::numeric_limits<double>::epsilon()) {
486 t_wk(
i) = t_w(
i) - t_k(
i);
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;
503 const double *v_ptr,
const int nb,
Range edges,
double *min_dist_ptr,
506 moab::Interface &moab(m_field.
get_moab());
511 auto get_point = [
i](
auto &t_w,
auto &t_delta,
auto t) {
513 t = std::max(0., std::min(1.,
t));
514 t_p(
i) = t_w(
i) +
t * t_delta(
i);
518 auto get_distance = [
i](
auto &t_p,
auto &t_n) {
520 t_dist_vector(
i) = t_p(
i) - t_n(
i);
521 return sqrt(t_dist_vector(
i) * t_dist_vector(
i));
524 for (
auto e : edges) {
528 CHKERR moab.get_connectivity(e, conn_fixed, num_nodes,
true);
530 CHKERR moab.get_coords(conn_fixed, num_nodes, &coords_fixed[0]);
532 &coords_fixed[0], &coords_fixed[1], &coords_fixed[2]);
534 &coords_fixed[3], &coords_fixed[4], &coords_fixed[5]);
537 t_edge_delta(
i) = t_f1(
i) - t_f0(
i);
540 v_ptr, v_ptr + 1, v_ptr + 2);
542 o_ptr, o_ptr + 1, o_ptr + 2);
547 colsest_segment_it = o_segments;
549 for (
int n = 0;
n != nb; ++
n) {
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) {
559 t_min_coords(
i) = t_p(
i);
561 *colsest_segment_it = e;
570 ++colsest_segment_it;
578 MatrixDouble &gauss_pts,
const int rule_ksi,
const int rule_eta) {
581 auto check_rule_edge = [](
int rule) {
585 "Wrong integration rule: %d", rule);
601 CHKERR check_rule_edge(rule_ksi);
602 CHKERR check_rule_edge(rule_eta);
606 gauss_pts.resize(3, nb_gauss_pts_ksi * nb_gauss_pts_eta,
false);
609 for (
size_t i = 0;
i != nb_gauss_pts_ksi; ++
i) {
612 for (
size_t j = 0;
j != nb_gauss_pts_eta; ++
j, ++gg) {
615 gauss_pts(0, gg) = ksi;
616 gauss_pts(1, gg) =
eta;
617 gauss_pts(2, gg) = wk;
620 if (gg != gauss_pts.size2())
627 MatrixDouble &gauss_pts,
const int rule_ksi,
const int rule_eta,
628 const int rule_zeta) {
631 auto check_rule_edge = [](
int rule) {
635 "Wrong integration rule: %d", rule);
651 CHKERR check_rule_edge(rule_ksi);
652 CHKERR check_rule_edge(rule_eta);
653 CHKERR check_rule_edge(rule_zeta);
658 gauss_pts.resize(4, nb_gauss_pts_ksi * nb_gauss_pts_eta * nb_gauss_pts_zeta,
662 for (
size_t i = 0;
i != nb_gauss_pts_ksi; ++
i) {
665 for (
size_t j = 0;
j != nb_gauss_pts_eta; ++
j) {
668 for (
size_t k = 0;
k != nb_gauss_pts_zeta; ++
k, ++gg) {
671 gauss_pts(0, gg) = ksi;
672 gauss_pts(1, gg) =
eta;
673 gauss_pts(2, gg) =
zeta;
674 gauss_pts(3, gg) = wk;
679 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.
double zeta
Viscous hardening.
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