v0.14.0
Classes | Public Types | Public Member Functions | Public Attributes | Private Member Functions | Private Attributes | List of all members
FractureMechanics::MWLSApprox Struct Reference

#include <users_modules/fracture_mechanics/src/MWLS.hpp>

Collaboration diagram for FractureMechanics::MWLSApprox:
[legend]

Classes

struct  OpMWLSBase
 
struct  OpMWLSCalculateBaseCoeffcientsAtGaussPtsTmpl
 
struct  OpMWLSMaterialStressLhs_Dx
 
struct  OpMWLSMaterialStressLhs_DX
 
struct  OpMWLSMaterialStressRhs
 Integrate force vector. More...
 
struct  OpMWLSRhoAtGaussPts
 Evaluate density at integration points. More...
 
struct  OpMWLSRhoAtGaussUsingPrecalulatedCoeffs
 
struct  OpMWLSRhoPostProcess
 
struct  OpMWLSSpatialStressLhs_DX
 
struct  OpMWLSSpatialStressRhs
 Integrate force vector. More...
 
struct  OpMWLSStressAndErrorsAtGaussPts
 Evaluate stress at integration points. More...
 
struct  OpMWLSStressAtGaussPts
 Evaluate stress at integration points. More...
 
struct  OpMWLSStressAtGaussUsingPrecalulatedCoeffs
 
struct  OpMWLSStressPostProcess
 

Public Types

typedef ublas::vector< double, ublas::bounded_array< double, 9 > > VecVals
 
typedef ublas::matrix< double, ublas::row_major, ublas::bounded_array< double, 27 > > MatDiffVals
 
typedef ublas::matrix< double, ublas::row_major, ublas::bounded_array< double, 81 > > MatDiffDiffVals
 
typedef OpMWLSCalculateBaseCoeffcientsAtGaussPtsTmpl< VolumeElementForcesAndSourcesCore > OpMWLSCalculateBaseCoeffcientsAtGaussPts
 
typedef OpMWLSCalculateBaseCoeffcientsAtGaussPtsTmpl< FaceElementForcesAndSourcesCoreOpMWLSCalculateBaseCoeffcientsForFaceAtGaussPts
 

Public Member Functions

 MWLSApprox (MoFEM::Interface &m_field, Vec F_lambda=PETSC_NULL, boost::shared_ptr< DofEntity > arc_length_dof=nullptr)
 
virtual ~MWLSApprox ()=default
 
bool getUseGlobalBaseAtMaterialReferenceConfiguration () const
 Get the Use Local Base At Material Reference Configuration object. More...
 
MoFEMErrorCode getMWLSOptions ()
 get options from command line for MWLS More...
 
MoFEMErrorCode loadMWLSMesh (std::string file_name)
 
MoFEMErrorCode getValuesToNodes (Tag th)
 get values form elements to nodes More...
 
MoFEMErrorCode getInfluenceNodes (double material_coords[3])
 Get node in influence domain. More...
 
MoFEMErrorCode calculateApproxCoeff (bool trim_influence_nodes, bool global_derivatives)
 Calculate approximation function coefficients. More...
 
MoFEMErrorCode evaluateApproxFun (double eval_point[3])
 evaluate approximation function at material point More...
 
MoFEMErrorCode evaluateFirstDiffApproxFun (double eval_point[3], bool global_derivatives)
 evaluate 1st derivative of approximation function at material point More...
 
MoFEMErrorCode evaluateSecondDiffApproxFun (double eval_point[3], bool global_derivatives)
 evaluate 2nd derivative of approximation function at material point More...
 
MoFEMErrorCode getTagData (Tag th)
 
MoFEMErrorCode getErrorAtNode (Tag th, double &total_stress_at_node, double &total_stress_error_at_node, double &hydro_static_error_at_node, double &deviatoric_error_at_node, double &total_energy, double &total_energy_error, const double young_modulus, const double poisson_ratio)
 Get the norm of the error at nod. More...
 
const VecValsgetDataApprox ()
 
const MatDiffValsgetDiffDataApprox ()
 
const MatDiffDiffValsgetDiffDiffDataApprox ()
 
MatrixDouble & getInvAB ()
 
std::vector< EntityHandle > & getInfluenceNodes ()
 
auto & getShortenDm ()
 

Public Attributes

MoFEM::InterfacemField
 
moab::Core mwlsCore
 
moab::Interface & mwlsMoab
 
Vec F_lambda
 
boost::weak_ptr< DofEntity > arcLengthDof
 
MatrixDouble stressAtGaussPts
 
MatrixDouble diffStressAtGaussPts
 
boost::shared_ptr< VectorDouble > rhoAtGaussPts
 
boost::shared_ptr< MatrixDouble > diffRhoAtGaussPts
 
boost::shared_ptr< MatrixDouble > diffDiffRhoAtGaussPts
 
MatrixDouble singularInitialDisplacement
 
MatrixDouble singularCurrentDisplacement
 
MatrixDouble mwlsMaterialPositions
 
MatrixDouble approxBasePoint
 
VectorDouble3 approxPointCoords
 
Range tetsInBlock
 
Range trisInBlock
 
Tag thMidNode
 
std::map< EntityHandle, std::vector< MatrixDouble > > invABMap
 Store Coefficient of base functions at integration points. More...
 
std::map< EntityHandle, std::vector< std::vector< EntityHandle > > > influenceNodesMap
 Influence nodes on elements on at integration points. More...
 
std::map< EntityHandle, std::vector< double > > dmNodesMap
 Store weight function radius at the nodes. More...
 
boost::function< VectorDouble(const double, const double, const double)> functionTestHack
 This is used to set up analytical function for testing approximation. More...
 

Private Member Functions

template<class T = FTensor::Tensor1<double, 3>>
MoFEMErrorCode calculateBase (const T &t_coords)
 
template<class T = FTensor::Tensor1<double, 3>>
MoFEMErrorCode calculateDiffBase (const T &t_coords, const double dm)
 
template<class T = FTensor::Tensor1<double, 3>>
MoFEMErrorCode calculateDiffDiffBase (const T &t_coords, const double dm)
 
double evalWeight (const double r)
 
double evalDiffWeight (const double r)
 
double evalDiffDiffWeight (const double r)
 

Private Attributes

boost::shared_ptr< BVHTree > myTree
 
int nbBasePolynomials
 
double dmFactor
 Relative value of weight function radius. More...
 
double shortenDm
 
int maxElemsInDMFactor
 
int maxThreeDepth
 Maximal three depths. More...
 
int splitsPerDir
 Split per direction in the tree. More...
 
PetscBool useGlobalBaseAtMaterialReferenceConfiguration
 
PetscBool useNodalData
 
PetscBool mwlsAnalyticalInternalStressTest
 
EntityHandle treeRoot
 
double maxEdgeL
 
Range levelNodes
 
Range levelEdges
 
Range levelTets
 
std::vector< EntityHandleinfluenceNodes
 
std::map< EntityHandle, std::array< double, 3 > > elemCentreMap
 
EntityHandle nearestInfluenceNode
 
ublas::symmetric_matrix< double, ublas::lower > A
 
ublas::symmetric_matrix< double, ublas::lower > testA
 
MatrixDouble invA
 
ublas::triangular_matrix< double, ublas::lower > L
 
MatrixDouble B
 
VectorDouble baseFun
 
VectorDouble approxFun
 
MatrixDouble invAB
 
VectorDouble baseFunInvA
 
ublas::symmetric_matrix< double, ublas::lower, ublas::row_major, DoubleAllocator > diffA [3]
 
ublas::symmetric_matrix< double, ublas::lower, ublas::row_major, DoubleAllocator > diffInvA [3]
 
MatrixDouble diffB [3]
 
VectorDouble diffBaseFun [3]
 
VectorDouble diffDiffBaseFun [9]
 
MatrixDouble diffApproxFun
 
MatrixDouble diffDiffApproxFun
 
MatrixDouble influenceNodesData
 
VecVals outData
 
double outHydroStaticError
 
VecVals outDeviatoricDifference
 
double outDeviatoricError
 
MatDiffVals outDiffData
 
MatDiffDiffVals outDiffDiffData
 
std::vector< EntityHandletreeTets
 
std::vector< EntityHandleleafsTetsVecLeaf
 
std::vector< doubleleafsTetsCentre
 
std::vector< EntityHandleleafsVec
 
std::vector< doubleleafsDist
 
std::vector< std::pair< double, EntityHandle > > distLeafs
 

Detailed Description

Moving least weighted square approximation (MWLS)

\[ u^h(x) = p_k(x)a_k(x)\\ J(x)=\frac{1}{2} \sum_i w(x-x_i)\left(p_j(x_i)a_j(x)-u_i\right)^2 = \\ J(x)=\frac{1}{2} \sum_i w(x-x_i)\left( p_j(x_i)a_j(x) p_k(x_i)a_k(x) - 2p_j(x_i)a_j(x)u_i + u_i^2 \right) \\ \frac{\partial J}{\partial a_j} = \sum_i w(x_i-x)\left(p_j(x_i)p_k(x_i) a_k(x) -p_j(x_j)u_i\right)\\ A_{jk}(x) = \sum_i w(x_i-x)p_j(x_i)p_k(x_i)\\ B_{ji}(x) = w(x-x_i) p_j(x_i)\\ A_{jk}(x) a_k(x) = B_{ji}(x)u_i\\ a_k(x) = \left( A^{-1}_{kj}(x) B_{ji}(x) \right) u_i\\ u^h(x) = p_k(x) \left( A^{-1}_{kj}(x) B_{ji}(x) \right) u_i \\ u^h(x) = \left[ p_k(x) \left( A^{-1}_{kj}(x) B_{ji}(x) \right) \right] u_i \\ \phi_i = p_k(x) \left( A^{-1}_{kj}(x) B_{ji}(x) \right) \\ \phi_{i,l} = p_{k,l}(x) \left( A^{-1}_{kj}(x) B_{ji}(x) \right)+ p_k(x) \left( A^{-1}_{kj,l}(x) B_{ji}(x) + A^{-1}_{kj}(x) B_{ji,l}(x) \right)\\ \phi_{i,l} = p_{k,l}(x) \left( A^{-1}_{kj}(x) B_{ji}(x) \right)+ p_k(x) \left( A^{-1}_{kj,l}(x) B_{ji}(x) + A^{-1}_{kj}(x) B_{ji,l}(x) \right)\\ \left(\mathbf{A}\mathbf{A}^{-1}\right)_{,l} = \mathbf{0} \to \mathbf{A}^{-1}_{,l} = -\mathbf{A}^{-1} \mathbf{A}_{,l} \mathbf{A}^{-1}\\ A_{jk,l} = \sum_i w_{,l}(x_i-x)p_j(x_i)p_k(x_i) \\ B_{ij,l} = w_{,l}(x_i-x)p_j(x_i) \]

Definition at line 53 of file MWLS.hpp.

Member Typedef Documentation

◆ MatDiffDiffVals

typedef ublas::matrix<double, ublas::row_major, ublas::bounded_array<double, 81> > FractureMechanics::MWLSApprox::MatDiffDiffVals

Definition at line 181 of file MWLS.hpp.

◆ MatDiffVals

typedef ublas::matrix<double, ublas::row_major, ublas::bounded_array<double, 27> > FractureMechanics::MWLSApprox::MatDiffVals

Definition at line 178 of file MWLS.hpp.

◆ OpMWLSCalculateBaseCoeffcientsAtGaussPts

Definition at line 296 of file MWLS.hpp.

◆ OpMWLSCalculateBaseCoeffcientsForFaceAtGaussPts

Definition at line 300 of file MWLS.hpp.

◆ VecVals

typedef ublas::vector<double, ublas::bounded_array<double, 9> > FractureMechanics::MWLSApprox::VecVals

Definition at line 175 of file MWLS.hpp.

Constructor & Destructor Documentation

◆ MWLSApprox()

FractureMechanics::MWLSApprox::MWLSApprox ( MoFEM::Interface m_field,
Vec  F_lambda = PETSC_NULL,
boost::shared_ptr< DofEntity arc_length_dof = nullptr 
)

Definition at line 36 of file MWLS.cpp.

38  : mField(m_field), mwlsMoab(mwlsCore), F_lambda(PETSC_NULL),
39  arcLengthDof(arc_length_dof), nbBasePolynomials(1), dmFactor(1),
42  useNodalData(PETSC_FALSE), mwlsAnalyticalInternalStressTest(PETSC_FALSE),
44 
45  if (!LogManager::checkIfChannelExist("MWLSWorld")) {
46  auto core_log = logging::core::get();
47 
48  core_log->add_sink(
49  LogManager::createSink(LogManager::getStrmWorld(), "MWLSWorld"));
50  core_log->add_sink(
51  LogManager::createSink(LogManager::getStrmSync(), "MWLSSync"));
52  core_log->add_sink(
53  LogManager::createSink(LogManager::getStrmSelf(), "MWLSSelf"));
54 
55  LogManager::setLog("MWLSWorld");
56  LogManager::setLog("MWLSSync");
57  LogManager::setLog("MWLSSelf");
58 
59  MOFEM_LOG_TAG("MWLSWorld", "MWLS");
60  MOFEM_LOG_TAG("MWLSSync", "MWLS");
61  MOFEM_LOG_TAG("MWLSSelf", "MWLS");
62  }
63 
64  MOFEM_LOG("MWLSWorld", Sev::noisy) << "MWLS created";
65 
66  rhoAtGaussPts = boost::make_shared<VectorDouble>();
67  diffRhoAtGaussPts = boost::make_shared<MatrixDouble>();
68  diffDiffRhoAtGaussPts = boost::make_shared<MatrixDouble>();
69 }

◆ ~MWLSApprox()

virtual FractureMechanics::MWLSApprox::~MWLSApprox ( )
virtualdefault

Member Function Documentation

◆ calculateApproxCoeff()

MoFEMErrorCode FractureMechanics::MWLSApprox::calculateApproxCoeff ( bool  trim_influence_nodes,
bool  global_derivatives 
)

Calculate approximation function coefficients.

Parameters
material_coordsmaterial position
Returns
error code

NOTE: This function remove form influenceNodes entities which are beyond influence radius.

Definition at line 576 of file MWLS.cpp.

577  {
579 
580  auto to_range = [&](auto ents) {
581  Range r;
582  r.insert_list(ents.begin(), ents.end());
583  return r;
584  };
585 
586  auto save_entities = [&](const std::string name, Range ents) {
588  if (!ents.empty()) {
589  EntityHandle meshset;
590  CHKERR mwlsMoab.create_meshset(MESHSET_SET, meshset);
591  CHKERR mwlsMoab.add_entities(meshset, ents);
592 
593  CHKERR mwlsMoab.write_file(name.c_str(), "VTK", "", &meshset, 1);
594  CHKERR mwlsMoab.delete_entities(&meshset, 1);
595  }
597  };
598 
599  // Determine points in influence volume
600  auto trim_nodes = [&](const double dm) {
603  FTensor::Tensor1<double, 3> t_approx_point(
605  const double dm2 = shortenDm * shortenDm;
606  VectorDouble node_coors_vec(3 * influenceNodes.size());
607  CHKERR mwlsMoab.get_coords(&*influenceNodes.begin(), influenceNodes.size(),
608  &*node_coors_vec.begin());
610  &node_coors_vec[0], &node_coors_vec[1], &node_coors_vec[2]);
611  decltype(influenceNodes) influence_nodes_tmp;
612  influence_nodes_tmp.reserve(influenceNodes.size());
613  for (auto node : influenceNodes) {
614  t_node_coords(i) -= t_approx_point(i);
615  const double r2 = t_node_coords(i) * t_node_coords(i);
616  if (r2 < dm2)
617  influence_nodes_tmp.emplace_back(node);
618  ++t_node_coords;
619  }
620  influenceNodes.swap(influence_nodes_tmp);
622  };
623 
624  auto eval_AB = [this, global_derivatives](const double dm) {
627 
628  A.resize(nbBasePolynomials, nbBasePolynomials, false);
629  B.resize(nbBasePolynomials, influenceNodes.size(), false);
630  A.clear();
631  B.clear();
632  for (int d = 0; d != 3; d++) {
633  diffA[d].resize(nbBasePolynomials, nbBasePolynomials, false);
634  diffB[d].resize(nbBasePolynomials, influenceNodes.size(), false);
635  diffA[d].clear();
636  diffB[d].clear();
637  }
638 
639  FTensor::Tensor1<double, 3> t_approx_point(
641 
642  VectorDouble nodes_coords(3 * influenceNodes.size());
643  CHKERR mwlsMoab.get_coords(&*influenceNodes.begin(), influenceNodes.size(),
644  &*nodes_coords.begin());
645  double *ptr = &*nodes_coords.begin();
647  ptr, ptr + 1, ptr + 2);
648 
649  double min_distance = -1;
651 
652  for (int ii = 0; ii != influenceNodes.size(); ++ii) {
653  t_node_coords(i) -= t_approx_point(i);
654  t_node_coords(i) /= shortenDm;
655  const double r = sqrt(t_node_coords(i) * t_node_coords(i));
656 
657  // Get the influence node that has minimum distance from the node in the
658  // projected mesh
659  if (!ii || r < min_distance) {
660  min_distance = r;
662  }
663 
665  if (r > 0) {
666  t_diff_r(i) = (1. / r) * t_node_coords(i) * (-1 / shortenDm);
667  } else {
668  t_diff_r(i) = 0;
669  }
670  // Weights
671  const double w = evalWeight(r);
672  calculateBase<FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>>(
673  t_node_coords);
674 
675  // iterate base functions
676  // A_kj = sum_ii w * b_k * b_j
677  //
678  // B_k = w * b_k
679  //
680  // diffA_kj = sum_ii diff_w * b_k * b_j + w * diff_b_k * b_j + w * b_k *
681  // diff_b_j
682  //
683  // diffB = diff_w * b_k + w * diff_bk
684  for (int k = 0; k != nbBasePolynomials; k++) {
685  const double wp = w * baseFun[k];
686  for (int j = 0; j <= k; j++) {
687  A(k, j) += wp * baseFun[j];
688  }
689  B(k, ii) = wp;
690  if (global_derivatives) {
691  const double diff_w_r = evalDiffWeight(r);
692  const double diff_wp_r = diff_w_r * baseFun[k];
693  for (int d = 0; d != 3; d++) {
694  const double diff_wp_r_dx = diff_wp_r * t_diff_r(d);
695  for (int j = 0; j <= k; j++) {
696  diffA[d](k, j) += diff_wp_r_dx * baseFun[j];
697  }
698  diffB[d](k, ii) = diff_wp_r_dx;
699  }
700  }
701  }
702  ++t_node_coords;
703  }
704  if (nearestInfluenceNode == 0) {
705  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
706  "Nearest node not found \n");
707  }
708 
710  };
711 
712  // Invert matrix A
713  auto invert_A = [this]() {
715 
716  auto solve = [&](const size_t nb_base_functions) {
717  A.resize(nb_base_functions, nb_base_functions, true);
718  L.resize(nb_base_functions, nb_base_functions, false);
719  if (cholesky_decompose(A, L)) {
720  return false;
721  } else {
722 
723  MatrixDouble inv_a(nb_base_functions, nb_base_functions, false);
724  for (int k = 0; k != nb_base_functions; k++) {
725  ublas::matrix_column<MatrixDouble> mr(inv_a, k);
726  mr(k) = 1;
727  cholesky_solve(L, mr, ublas::lower());
728  }
729 
730  if (nb_base_functions == nbBasePolynomials) {
731  invA.swap(inv_a);
732  } else {
733  invA.resize(nbBasePolynomials, nbBasePolynomials, false);
734  invA.clear();
735  for (size_t i = 0; i != nb_base_functions; ++i)
736  for (size_t j = 0; j != nb_base_functions; ++j)
737  invA(i, j) = inv_a(i, j);
738  }
739 
740  return true;
741  }
742  };
743 
744  auto throw_error = [&]() {
746  FTensor::Tensor1<double, 3> t_approx_point(
748  cerr << "Point: " << t_approx_point << endl;
749  cerr << "Matrix A: " << A << endl;
751  "Failed to invert matrix - solution could be "
752  "to increase dmFactor to bigger value, e.g. -mwls_dm 4, "
753  "or reduce of number of base fuctions "
754  "-mwls_number_of_base_functions 1. The number of nodes "
755  "found in "
756  "radius is %u",
757  influenceNodes.size());
759  };
760 
761  auto solve_one = [&]() {
763  if (!solve(1))
764  CHKERR throw_error();
766  };
767 
768  testA = A;
769  if (!solve(nbBasePolynomials)) {
770  if (nbBasePolynomials == 10) {
771  if (!solve(4))
772  CHKERR solve_one();
773  } else if (nbBasePolynomials == 4) {
774  CHKERR solve_one();
775  } else {
776  CHKERR throw_error();
777  }
778  }
779 
781  };
782 
783  auto calculate_shape_functions_coefficients = [this]() {
785  // Calculate base functions (shape functions) coefficients
786  invAB.resize(invA.size1(), B.size2(), false);
787  noalias(invAB) = prod(invA, B);
789  };
790 
791  // CHKERR trim_nodes(shortenDm);
792  // CHKERR save_entities("tree_tets_trim.vtk", to_range(influenceNodes));
793 
794  CHKERR eval_AB(shortenDm);
795  CHKERR invert_A();
796  CHKERR calculate_shape_functions_coefficients();
797 
799 }

◆ calculateBase()

template<class T = FTensor::Tensor1<double, 3>>
MoFEMErrorCode FractureMechanics::MWLSApprox::calculateBase ( const T &  t_coords)
inlineprivate

Definition at line 675 of file MWLS.hpp.

675  {
677 
678  baseFun.resize(nbBasePolynomials, false);
679  baseFun.clear();
680 
681  baseFun[0] = 1;
682  if (nbBasePolynomials > 1) {
683  baseFun[1] = t_coords(0); // x
684  baseFun[2] = t_coords(1); // y
685  baseFun[3] = t_coords(2); // z
686  }
687 
688  if (nbBasePolynomials > 4) {
689  baseFun[4] = t_coords(0) * t_coords(0); // x^2
690  baseFun[5] = t_coords(1) * t_coords(1); // y^2
691  baseFun[6] = t_coords(2) * t_coords(2); // z^2
692  baseFun[7] = t_coords(0) * t_coords(1); // xy
693  baseFun[8] = t_coords(0) * t_coords(2); // xz
694  baseFun[9] = t_coords(1) * t_coords(2); // yz
695  }
696 
698  }

◆ calculateDiffBase()

template<class T = FTensor::Tensor1<double, 3>>
MoFEMErrorCode FractureMechanics::MWLSApprox::calculateDiffBase ( const T &  t_coords,
const double  dm 
)
inlineprivate

Definition at line 701 of file MWLS.hpp.

701  {
703 
704  // Local base functions
705  for (int d = 0; d != 3; ++d) {
706  diffBaseFun[d].resize(nbBasePolynomials, false);
707  diffBaseFun[d].clear();
708  }
709 
710  for (int d = 0; d != 3; d++) {
711  diffBaseFun[d][0] = 0;
712  }
713 
714  if (nbBasePolynomials > 1) {
715  // derivative dx
716  diffBaseFun[0][1] = 1. / dm;
717  // derivative dy
718  diffBaseFun[1][2] = 1. / dm;
719  // derivative dz
720  diffBaseFun[2][3] = 1. / dm;
721  }
722 
723  if (nbBasePolynomials > 4) {
724  // dx derivative
725  diffBaseFun[0][4] = 2 * t_coords(0) / dm;
726  diffBaseFun[0][5] = 0;
727  diffBaseFun[0][6] = 0;
728  diffBaseFun[0][7] = t_coords(1) / dm;
729  diffBaseFun[0][8] = t_coords(2) / dm;
730  diffBaseFun[0][9] = 0;
731  // dy derivative
732  diffBaseFun[1][4] = 0;
733  diffBaseFun[1][5] = 2 * t_coords(1) / dm;
734  diffBaseFun[1][6] = 0;
735  diffBaseFun[1][7] = t_coords(0) / dm;
736  diffBaseFun[1][8] = 0;
737  diffBaseFun[1][9] = t_coords(2) / dm;
738  // dz derivative
739  diffBaseFun[2][4] = 0;
740  diffBaseFun[2][5] = 0;
741  diffBaseFun[2][6] = 2 * t_coords(2) / dm;
742  diffBaseFun[2][7] = 0;
743  diffBaseFun[2][8] = t_coords(0) / dm;
744  diffBaseFun[2][9] = t_coords(1) / dm;
745  }
747  }

◆ calculateDiffDiffBase()

template<class T = FTensor::Tensor1<double, 3>>
MoFEMErrorCode FractureMechanics::MWLSApprox::calculateDiffDiffBase ( const T &  t_coords,
const double  dm 
)
inlineprivate

Definition at line 750 of file MWLS.hpp.

750  {
752 
753  // Local base functions
754  for (int d = 0; d != 9; ++d) {
755  diffDiffBaseFun[d].resize(nbBasePolynomials, false);
756  diffDiffBaseFun[d].clear();
757  }
758 
759  // if (nbBasePolynomials > 1) {
760  // }
761 
762  if (nbBasePolynomials > 4) {
763  const double dm2 = dm * dm;
764  // set the non-zero values
765  // dxd? derivative
766  diffDiffBaseFun[0][4] = 2. / dm2; // dxdx
767  diffDiffBaseFun[1][7] = 1. / dm2; // dxdy
768  diffDiffBaseFun[2][8] = 1. / dm2; // dxdz
769  // dyd? derivative
770  diffDiffBaseFun[4][5] = 2. / dm2; // dydy
771  diffDiffBaseFun[3][7] = 1. / dm2; // dydx
772  diffDiffBaseFun[5][9] = 1. / dm2; // dydz
773  // dzd? derivative
774  diffDiffBaseFun[8][6] = 2. / dm2; // dzdz
775  diffDiffBaseFun[6][8] = 1. / dm2; // dzdx
776  diffDiffBaseFun[7][9] = 1. / dm2; // dzdy
777  }
779  }

◆ evalDiffDiffWeight()

double FractureMechanics::MWLSApprox::evalDiffDiffWeight ( const double  r)
inlineprivate

Definition at line 812 of file MWLS.hpp.

812  {
813  const double rr = r * r;
814  if (r < 1) {
815  return -36 * rr + 48 * r - 12;
816  } else {
817  return 0;
818  }
819  }

◆ evalDiffWeight()

double FractureMechanics::MWLSApprox::evalDiffWeight ( const double  r)
inlineprivate

Definition at line 797 of file MWLS.hpp.

797  {
798  const double rr = r * r;
799  const double rrr = rr * r;
800  if (r < 1) {
801  return -12 * r + 24 * rr - 12 * rrr;
802  } else {
803  return 0;
804  }
805  // if(r<=0.5) {
806  // return -8*r+12*rr;
807  // } else if(r<=1) {
808  // return -4+8*r-4*rr;
809  // } else return 0;
810  }

◆ evaluateApproxFun()

MoFEMErrorCode FractureMechanics::MWLSApprox::evaluateApproxFun ( double  eval_point[3])

evaluate approximation function at material point

Parameters
material_coordsmaterial position
Returns
error code

Definition at line 801 of file MWLS.cpp.

801  {
803  // Determine points in influence volume
804 
805  auto cal_base_functions_at_eval_point = [this, eval_point](const double dm) {
807 
808  FTensor::Tensor1<double, 3> t_approx_point(
810  FTensor::Tensor1<double, 3> t_eval_point(eval_point[0], eval_point[1],
811  eval_point[2]);
812 
815  t_delta(i) = t_eval_point(i) - t_approx_point(i);
816  t_delta(i) /= dm;
817 
818  CHKERR calculateBase(t_delta);
820  };
821 
822  auto calculate_shape_functions = [this]() {
824  // Calculate approximation functions at eval point
825  approxFun.resize(influenceNodes.size(), false);
826  noalias(approxFun) = prod(baseFun, invAB);
828  };
829 
830  CHKERR cal_base_functions_at_eval_point(shortenDm);
831  CHKERR calculate_shape_functions();
832 
834 }

◆ evaluateFirstDiffApproxFun()

MoFEMErrorCode FractureMechanics::MWLSApprox::evaluateFirstDiffApproxFun ( double  eval_point[3],
bool  global_derivatives 
)

evaluate 1st derivative of approximation function at material point

Parameters
global_derivativesdecide whether global derivative will be calculated
Returns
error code

Definition at line 836 of file MWLS.cpp.

837  {
839 
840  auto cal_diff_base_functions_at_eval_point = [this,
841  eval_point](const double dm) {
843 
844  FTensor::Tensor1<double, 3> t_approx_point(
846  FTensor::Tensor1<double, 3> t_eval_point(eval_point[0], eval_point[1],
847  eval_point[2]);
848 
851  t_delta(i) = t_eval_point(i) - t_approx_point(i);
852  t_delta(i) /= dm;
853 
854  CHKERR calculateDiffBase(t_delta, dm);
856  };
857 
858  auto calculate_global_shape_functions_derivatives = [this]() {
860 
861  // Calculate derivatives of base functions
862  VectorDouble tmp;
863  tmp.resize(nbBasePolynomials, false);
864  baseFunInvA.resize(invA.size2(), false);
865  noalias(baseFunInvA) = prod(baseFun, invA);
866 
867  for (int d = 0; d != 3; d++) {
869  diffApproxFun.resize(3, influenceNodes.size(), false);
870  }
871 
872  for (int d = 0; d != 3; d++) {
873  MatrixDouble a = prod(diffA[d], invA);
874  noalias(diffInvA[d]) = -prod(invA, a);
875  }
876 
877  // Calculate line derivatives for each coordinate direction
878  for (int d = 0; d != 3; ++d) {
879 
880  ublas::matrix_row<MatrixDouble> mr(diffApproxFun, d);
881 
882  // b,x * A^-1 B
883  noalias(mr) = prod(diffBaseFun[d], invAB);
884 
885  // b * (A^-1 A,x A^-1) * B
886  noalias(tmp) = prod(baseFun, diffInvA[d]);
887  mr += prod(tmp, B);
888 
889  // b * A^-1 * B,x
890  mr += prod(baseFunInvA, diffB[d]);
891  }
892 
894  };
895 
896  auto calculate_local_shape_functions_derivatives = [this]() {
898 
899  // Calculate derivatives of base functions
900  diffApproxFun.resize(3, influenceNodes.size(), false);
901  // Calculate line derivatives for each coordinate direction
902  for (int d = 0; d != 3; ++d) {
903  ublas::matrix_row<MatrixDouble> mr(diffApproxFun, d);
904  // b,x * A^-1 B
905  noalias(mr) = prod(diffBaseFun[d], invAB);
906  }
907 
909  };
910 
911  CHKERR cal_diff_base_functions_at_eval_point(shortenDm);
912  if (global_derivatives)
913  CHKERR calculate_global_shape_functions_derivatives();
914  else
915  CHKERR calculate_local_shape_functions_derivatives();
916 
918 }

◆ evaluateSecondDiffApproxFun()

MoFEMErrorCode FractureMechanics::MWLSApprox::evaluateSecondDiffApproxFun ( double  eval_point[3],
bool  global_derivatives 
)

evaluate 2nd derivative of approximation function at material point

Parameters
global_derivativesdecide whether global derivative will be calculated
Returns
error code

Definition at line 921 of file MWLS.cpp.

922  {
924 
925  auto cal_diff_diff_base_functions_at_eval_point =
926  [this, eval_point](const double dm) {
928 
929  FTensor::Tensor1<double, 3> t_approx_point(
931  FTensor::Tensor1<double, 3> t_eval_point(eval_point[0], eval_point[1],
932  eval_point[2]);
933 
936  t_delta(i) = t_eval_point(i) - t_approx_point(i);
937  t_delta(i) /= dm;
938 
939  CHKERR calculateDiffDiffBase(t_delta, dm);
941  };
942 
943  auto calculate_global_shape_functions_second_derivatives = [this]() {
945  // to be implemented
947  };
948 
949  auto calculate_local_shape_functions_second_derivatives = [this]() {
951 
952  // Calculate derivatives of base functions
953  // for (int d = 0; d != 3; d++) {
954  diffDiffApproxFun.resize(9, influenceNodes.size(), false);
955  // }
956  // Calculate line derivatives for each coordinate direction
957  for (int n = 0; n != 3; ++n) {
958  for (int d = 0; d != 3; ++d) {
959  const int indx = d + 3 * n;
960  ublas::matrix_row<MatrixDouble> mr(diffDiffApproxFun, indx);
961  // b,xy * A^-1 B
962  noalias(mr) = prod(diffDiffBaseFun[indx], invAB);
963  }
964  }
966  };
967 
968  CHKERR cal_diff_diff_base_functions_at_eval_point(shortenDm);
969  if (global_derivatives) {
970  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented yet");
971  // CHKERR calculate_global_shape_functions_second_derivatives();
972  } else
973  CHKERR calculate_local_shape_functions_second_derivatives();
974 
976 }

◆ evalWeight()

double FractureMechanics::MWLSApprox::evalWeight ( const double  r)
inlineprivate

Definition at line 781 of file MWLS.hpp.

781  {
782  const double rr = r * r;
783  const double rrr = rr * r;
784  const double r4 = rrr * r;
785  if (r < 1) {
786  return 1 - 6 * rr + 8 * rrr - 3 * r4;
787  } else {
788  return 0;
789  }
790  // if(r<=0.5) {
791  // return 2./3-4*rr+4*rrr;
792  // } else if(r<=1) {
793  // return 4./3-4*r+4*rr-(4./3)*rrr;
794  // } else return 0;
795  }

◆ getDataApprox()

const MWLSApprox::VecVals & FractureMechanics::MWLSApprox::getDataApprox ( )

Get values at point

Returns
error code

Definition at line 992 of file MWLS.cpp.

992  {
993  noalias(outData) = prod(approxFun, influenceNodesData);
994  return outData;
995 }

◆ getDiffDataApprox()

const MWLSApprox::MatDiffVals & FractureMechanics::MWLSApprox::getDiffDataApprox ( )

Get derivatives of values at points

Returns
error code

Definition at line 997 of file MWLS.cpp.

997  {
998  noalias(outDiffData) = prod(diffApproxFun, influenceNodesData);
999  return outDiffData;
1000 }

◆ getDiffDiffDataApprox()

const MWLSApprox::MatDiffDiffVals & FractureMechanics::MWLSApprox::getDiffDiffDataApprox ( )

Get second derivatives of values at points

Returns
error code

Definition at line 1002 of file MWLS.cpp.

1002  {
1004  return outDiffDiffData;
1005 }

◆ getErrorAtNode()

MoFEMErrorCode FractureMechanics::MWLSApprox::getErrorAtNode ( Tag  th,
double total_stress_at_node,
double total_stress_error_at_node,
double hydro_static_error_at_node,
double deviatoric_error_at_node,
double total_energy,
double total_energy_error,
const double  young_modulus,
const double  poisson_ratio 
)

Get the norm of the error at nod.

Parameters
th
total_stress_at_node
total_stress_error_at_node
hydro_static_error_at_node
deviatoric_error_at_node
total_energy
total_energy_error
young_modulus
poisson_ratio
Returns
MoFEMErrorCode

Definition at line 1007 of file MWLS.cpp.

1011  {
1013 
1014  int length;
1015  CHKERR mwlsMoab.tag_get_length(th, length);
1016 
1021 
1022  auto get_compliance_matrix = [&]() {
1024  t_C(i, j, k, l) = 0;
1025 
1026  const double c = poisson_ratio / young_modulus;
1027  t_C(0, 0, 0, 0) = 1. / young_modulus;
1028  t_C(0, 0, 1, 1) = -c;
1029  t_C(0, 0, 2, 2) = -c;
1030 
1031  t_C(1, 1, 0, 0) = -c;
1032  t_C(1, 1, 1, 1) = 1. / young_modulus;
1033  t_C(1, 1, 2, 2) = -c;
1034 
1035  t_C(2, 2, 0, 0) = -c;
1036  t_C(2, 2, 1, 1) = -c;
1037  t_C(2, 2, 2, 2) = 1. / young_modulus;
1038 
1039  const double d = (1 + poisson_ratio) / young_modulus;
1040  t_C(0, 1, 0, 1) = d;
1041  t_C(0, 2, 0, 2) = d;
1042  t_C(1, 2, 1, 2) = d;
1043 
1044  return t_C;
1045  };
1046 
1047  auto t_C = get_compliance_matrix();
1048 
1049  if (length != outData.size())
1050  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1051  "Tag length is not consistent with outData size. \nMost probably "
1052  "getTagData has not been invoked yet!\n");
1053 
1054  VectorDouble val_at_nearest_influence_node;
1055  val_at_nearest_influence_node.resize(length, false);
1056  val_at_nearest_influence_node.clear();
1057  CHKERR mwlsMoab.tag_get_data(th, &nearestInfluenceNode, 1,
1058  &*val_at_nearest_influence_node.begin());
1059 
1060  // SIXX 0 (00) SIYY 1 (11) SIZZ 2 (22) SIXY 3 (01) SIXZ 4 (02) SIYZ 5 (12)
1061  // T* d00, T* d01, T* d02, T* d11, T* d12, T* d22
1062  auto get_symm_tensor = [](auto &m) {
1064 
1065  &m(0), &m(3), &m(4),
1066 
1067  &m(1), &m(5),
1068 
1069  &m(2)
1070 
1071  );
1072  return t;
1073  };
1074 
1075  auto t_total_stress_error = FTensor::Tensor2_symmetric<double, 3>();
1076  auto t_stress_at_nearset_node =
1077  get_symm_tensor(val_at_nearest_influence_node);
1078  auto t_mwls_stress = get_symm_tensor(outData);
1079 
1080  auto get_energy_norm = [&](auto &t_s) {
1081  return 0.5 * t_s(i, j) * (t_C(i, j, k, l) * t_s(k, l));
1082  };
1083  auto get_stress_nrm2 = [i, j](auto &t_s) { return t_s(i, j) * t_s(i, j); };
1084 
1085  total_stress_at_node = get_stress_nrm2(t_stress_at_nearset_node);
1086  total_energy = get_energy_norm(t_stress_at_nearset_node);
1087 
1088  t_total_stress_error(i, j) =
1089  t_mwls_stress(i, j) - t_stress_at_nearset_node(i, j);
1090  total_stress_error_at_node = get_stress_nrm2(t_total_stress_error);
1091  total_energy_error = get_energy_norm(t_total_stress_error);
1092 
1093  constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
1094  // calculate error of the hydro static part of the stress
1095  double hydro_static_difference =
1096  (t_kd(i, j) * t_total_stress_error(i, j)) / 3.;
1097  hydro_static_error_at_node =
1098  hydro_static_difference * hydro_static_difference;
1099 
1100  // calculate error of the deviator part of the stress
1101  auto t_dev_error = FTensor::Tensor2_symmetric<double, 3>();
1102  t_dev_error(i, j) =
1103  t_total_stress_error(i, j) - hydro_static_difference * t_kd(i, j);
1104  deviatoric_error_at_node = get_stress_nrm2(t_dev_error);
1105 
1107 }

◆ getInfluenceNodes() [1/2]

std::vector<EntityHandle>& FractureMechanics::MWLSApprox::getInfluenceNodes ( )
inline

Definition at line 605 of file MWLS.hpp.

605  {
606  return influenceNodes;
607  }

◆ getInfluenceNodes() [2/2]

MoFEMErrorCode FractureMechanics::MWLSApprox::getInfluenceNodes ( double  material_coords[3])

Get node in influence domain.

Parameters
material_coordsmaterial position
Returns
error code

Definition at line 373 of file MWLS.cpp.

373  {
374  constexpr double eps_overlap = 0.01;
375  const size_t max_nb_elems = maxElemsInDMFactor * nbBasePolynomials + 1;
377 
378  using DistEntPair = std::pair<double, EntityHandle>;
379  using DistEntMap = multi_index_container<
380  DistEntPair,
381  indexed_by<
382 
383  ordered_non_unique<member<DistEntPair, double, &DistEntPair::first>>
384 
385  >>;
386  DistEntMap dist_map;
387 
388  // FIXME: That is potentisl source of problems, should be sparate
389  // function for setting approximation base point, or simply should be set
390  // outside of that function. Other functions will behave depending how
391  // approxPointCoords is set. And here is somehow hidden from the view,
392  // i.e. easy to overlook this. TODO: create a new function like:
393  // MoFEMErrorCode setApproxBasePt(double material_coords[3]) ?
394 
395  for (int dd : {0, 1, 2})
396  approxPointCoords[dd] = material_coords[dd];
397 
398  if (myTree) {
399  shortenDm = maxEdgeL * dmFactor; // Influence radius
400 
401  auto find_leafs = [&](const auto dm) {
402  leafsVec.clear();
403  leafsDist.clear();
404  CHKERR myTree->distance_search(material_coords, shortenDm, leafsVec,
405  1.0e-10, 1.0e-6, &leafsDist);
406  distLeafs.clear();
407  distLeafs.reserve(leafsDist.size());
408  for (size_t i = 0; i != leafsVec.size(); ++i) {
409  distLeafs.emplace_back(leafsDist[i], leafsVec[i]);
410  }
411  std::sort(
412  distLeafs.begin(), distLeafs.end(),
413  [](const auto &a, const auto &b) { return a.first < b.first; });
414  return distLeafs;
415  };
416 
417  auto get_influence_nodes = [&](const auto &tets) {
418  std::vector<EntityHandle> influence_nodes;
419  if (!useNodalData) {
420  influence_nodes.resize(tets.size());
421  CHKERR mwlsMoab.tag_get_data(thMidNode, &*tets.begin(), tets.size(),
422  &*influence_nodes.begin());
423  } else {
424  Range nodes;
425  CHKERR mwlsMoab.get_connectivity(&*tets.begin(), tets.size(), nodes,
426  true);
427  influence_nodes.clear();
428  influence_nodes.reserve(nodes.size());
429  for (auto n : nodes)
430  influence_nodes.emplace_back(n);
431  }
432  return influence_nodes;
433  };
434 
435  auto to_range = [&](auto ents) {
436  Range r;
437  r.insert_list(ents.begin(), ents.end());
438  return r;
439  };
440 
441  auto save_entities = [&](const std::string name, Range ents) {
443  if (!ents.empty()) {
444  EntityHandle meshset;
445  CHKERR mwlsMoab.create_meshset(MESHSET_SET, meshset);
446  CHKERR mwlsMoab.add_entities(meshset, ents);
447 
448  CHKERR mwlsMoab.write_file(name.c_str(), "VTK", "", &meshset, 1);
449  CHKERR mwlsMoab.delete_entities(&meshset, 1);
450  }
452  };
453 
454  auto leafs = find_leafs(shortenDm);
455 
456  treeTets.clear();
457  std::vector<EntityHandle> leafs_tets;
458  moab::BoundBox box;
459 
460  auto get_dm2_from_influence_points = [&]() {
462  VectorDouble node_coors_vec(3 * influenceNodes.size());
463  CHKERR mwlsMoab.get_coords(&*influenceNodes.begin(),
464  influenceNodes.size(),
465  &*node_coors_vec.begin());
466  FTensor::Tensor1<double, 3> t_approx_point(
467  material_coords[0], material_coords[1], material_coords[2]);
469  &node_coors_vec[0], &node_coors_vec[1], &node_coors_vec[2]);
470  double dm2 = 0;
471  for (auto node : influenceNodes) {
472  t_node_coords(i) -= t_approx_point(i);
473  const double r2 = t_node_coords(i) * t_node_coords(i);
474  if (dm2 < r2)
475  dm2 = r2;
476  ++t_node_coords;
477  }
478  return dm2;
479  };
480 
481  const double shorten_dm2 = shortenDm * shortenDm;
482  double set_dm2 = shorten_dm2;
483  auto min_max_dist_it = dist_map.get<0>().begin();
485 
486  for (auto l : leafs) {
487  if (dist_map.size() < max_nb_elems ||
488  ((l.first * l.first) <
489  (set_dm2 + std::numeric_limits<double>::epsilon()))) {
490 
491  leafsTetsVecLeaf.clear();
492  CHKERR mwlsMoab.get_entities_by_dimension(l.second, 3,
493  leafsTetsVecLeaf, false);
494 
495  leafsTetsCentre.resize(3 * leafsTetsVecLeaf.size());
496  {
499  &leafsTetsCentre[2 * leafsTetsVecLeaf.size()]};
500 
501  for (auto tet : leafsTetsVecLeaf) {
502  const auto &c = elemCentreMap[tet];
503  for (auto ii : {0, 1, 2}) {
504  t_c(ii) = c[ii] - material_coords[ii];
505  }
506  ++t_c;
507  }
508  }
509 
510  {
513  &leafsTetsCentre[2 * leafsTetsVecLeaf.size()]};
514 
515  for (auto tet : leafsTetsVecLeaf) {
516 
517  if (dist_map.size() > max_nb_elems)
518  set_dm2 = std::min(min_max_dist_it->first, shorten_dm2);
519 
520  if (t_c(0) < set_dm2 && t_c(1) < set_dm2 && t_c(2) < set_dm2) {
521  double dm2 = t_c(i) * t_c(i);
522  if (dm2 < set_dm2) {
523  min_max_dist_it =
524  dist_map.get<0>().emplace(std::make_pair(dm2, tet)).first;
525  if (dist_map.size() > max_nb_elems) {
526  int nb = std::distance(dist_map.get<0>().begin(),
527  min_max_dist_it);
528  for (; nb < max_nb_elems; ++nb)
529  ++min_max_dist_it;
530  }
531  }
532  }
533 
534  ++t_c;
535  }
536  }
537  } else {
538  break;
539  }
540  }
541 
542  treeTets.clear();
543  treeTets.reserve(dist_map.size());
544  const double dm2 = set_dm2 + std::numeric_limits<double>::epsilon();
545  auto hi_it = dist_map.get<0>().upper_bound(dm2);
546  for (auto it = dist_map.get<0>().begin(); it != hi_it; ++it)
547  treeTets.emplace_back(it->second);
548 
549  if (treeTets.empty()) {
550  for (auto &it : dist_map)
551  MOFEM_LOG("MWLSWorld", Sev::error)
552  << "dm map " << it.first << " " << it.second;
553 
554  MOFEM_LOG("MWLSWorld", Sev::error) << "leafs found " << leafs.size();
555  MOFEM_LOG("MWLSWorld", Sev::error)
556  << "dist map size " << dist_map.size();
557  MOFEM_LOG("MWLSWorld", Sev::error) << "dm2 " << dm2;
558  MOFEM_LOG("MWLSWorld", Sev::error) << "shorten dm is " << shortenDm;
559  MOFEM_LOG_C("MWLSWorld", Sev::error, "point: ( %g %g %g )",
561  approxPointCoords(2));
562  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No leafs found.");
563  }
564 
565  influenceNodes = get_influence_nodes(treeTets);
566  const auto dm2_from_influence_pts = get_dm2_from_influence_points();
567  shortenDm = (1 + eps_overlap) * sqrt(dm2_from_influence_pts);
568 
569  // CHKERR save_entities("tree_tets.vtk", to_range(influenceNodes));
570  } else {
571  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "kd-three not build");
572  }
574 }

◆ getInvAB()

MatrixDouble& FractureMechanics::MWLSApprox::getInvAB ( )
inline

Definition at line 604 of file MWLS.hpp.

604 { return invAB; }

◆ getMWLSOptions()

MoFEMErrorCode FractureMechanics::MWLSApprox::getMWLSOptions ( )

get options from command line for MWLS

Returns
error code

Definition at line 71 of file MWLS.cpp.

71  {
73  ierr = PetscOptionsBegin(mField.get_comm(), "",
74  "Get options for moving least square approximation",
75  "none");
76  CHKERRQ(ierr);
77  CHKERR PetscOptionsScalar("-mwls_dm",
78  "edge length scaling for influence radius", "",
79  dmFactor, &dmFactor, PETSC_NULL);
80  MOFEM_LOG_C("MWLSWorld", Sev::inform, "### Input parameter: -mwls_dm %6.4e",
81  dmFactor);
82 
83  CHKERR PetscOptionsInt(
84  "-mwls_number_of_base_functions",
85  "1 = constant, 4 = linear, 10 = quadratic approximation", "",
86  nbBasePolynomials, &nbBasePolynomials, PETSC_NULL);
87  MOFEM_LOG_C("MWLSWorld", Sev::inform,
88  "### Input parameter: -mwls_number_of_base_functions %d",
90 
91  CHKERR PetscOptionsInt("-mwls_max_elems_factor",
92  "Set max elements factor max_nb_elems = "
93  "maxElemsInDMFactor * nbBasePolynomials",
95  PETSC_NULL);
96  MOFEM_LOG_C("MWLSWorld", Sev::inform,
97  "### Input parameter: -mwls_max_elems_factor %d",
99 
100  CHKERR PetscOptionsBool(
101  "-mwls_use_global_base_at_reference_configuration",
102  "true local mwls base functions at reference configuration are used", "",
105  CHKERRQ(ierr);
106 
107  CHKERR PetscOptionsBool(
108  "-mwls_use_nodal_data",
109  "true nodal data values are used (without averaging to the midpoint)", "",
110  useNodalData, &useNodalData, NULL);
111  MOFEM_LOG_C("MWLSWorld", Sev::verbose,
112  "### Input parameter: -mwls_use_nodal_data %d", useNodalData);
113 
114  CHKERR PetscOptionsBool(
115  "-mwls_analytical_internal_stress",
116  "use analytical strain field for internal stress mwls test", "",
118  NULL);
119  MOFEM_LOG_C("MWLSWorld", Sev::verbose,
120  "### Input parameter: -test_mwls_internal_stress_analytical %d",
122 
123  maxThreeDepth = 30;
124  CHKERR PetscOptionsInt("-mwls_tree_depth", "maximal three depths", "",
125  maxThreeDepth, &maxThreeDepth, PETSC_NULL);
126  MOFEM_LOG("MWLSWorld", Sev::verbose)
127  << "Maximal MWLS three depths " << maxThreeDepth;
128  splitsPerDir = 12;
129  CHKERR PetscOptionsInt("-mwls_splits_per_dir", "splits per direction", "",
130  splitsPerDir, &splitsPerDir, PETSC_NULL);
131  MOFEM_LOG("MWLSWorld", Sev::verbose)
132  << "Splits per direction " << splitsPerDir;
133 
134  ierr = PetscOptionsEnd();
135  CHKERRQ(ierr);
137 }

◆ getShortenDm()

auto& FractureMechanics::MWLSApprox::getShortenDm ( )
inline

Definition at line 608 of file MWLS.hpp.

608 { return shortenDm; }

◆ getTagData()

MoFEMErrorCode FractureMechanics::MWLSApprox::getTagData ( Tag  th)

Get tag data on influence nodes

Parameters
thtag
Returns
error code

Definition at line 978 of file MWLS.cpp.

978  {
980  int length;
981  CHKERR mwlsMoab.tag_get_length(th, length);
982  influenceNodesData.resize(influenceNodes.size(), length, false);
983  CHKERR mwlsMoab.tag_get_data(th, &*influenceNodes.begin(),
984  influenceNodes.size(),
985  &influenceNodesData(0, 0));
986  outData.resize(length, false);
987  outDiffData.resize(3, length, false);
988  outDiffDiffData.resize(9, length, false);
990 }

◆ getUseGlobalBaseAtMaterialReferenceConfiguration()

bool FractureMechanics::MWLSApprox::getUseGlobalBaseAtMaterialReferenceConfiguration ( ) const
inline

Get the Use Local Base At Material Reference Configuration object.

If true base is calculated at initail integration point coordinates and base function evaluated at another point. Local derivatives are uses, since weight function is at initial point.

Returns
true
false

Definition at line 76 of file MWLS.hpp.

76  {
77  return (useGlobalBaseAtMaterialReferenceConfiguration == PETSC_TRUE);
78  }

◆ getValuesToNodes()

MoFEMErrorCode FractureMechanics::MWLSApprox::getValuesToNodes ( Tag  th)

get values form elements to nodes

Parameters
thtag with approximated data
thtag with approximated data
Returns
error code

Definition at line 275 of file MWLS.cpp.

275  {
277  if (useNodalData)
279 
280  CHKERR mwlsMoab.tag_get_handle("MID_NODE", 1, MB_TYPE_HANDLE, thMidNode,
281  MB_TAG_CREAT | MB_TAG_DENSE);
282 
283  int length;
284  CHKERR mwlsMoab.tag_get_length(th, length);
285 
286  VectorDouble vals(length);
287  vals.clear();
288 
289  if (mwlsAnalyticalInternalStressTest == PETSC_FALSE) {
290 
291  string name;
292  CHKERR mwlsMoab.tag_get_name(th, name);
293 
294  std::function<MoFEMErrorCode(EntityHandle)> get_data_to_node;
295  if (name.compare("RHO") == 0) {
296  get_data_to_node = [this, th, &vals](EntityHandle tet) {
298  Range nodes;
299  double avg_val = 0;
300  CHKERR mwlsMoab.get_connectivity(&tet, 1, nodes, true);
301  for (auto &n : nodes) {
302  double tag_data;
303  CHKERR mwlsMoab.tag_get_data(th, &n, 1, &tag_data);
304  avg_val += tag_data;
305  }
306  avg_val /= nodes.size();
307  vals[0] = avg_val;
309  };
310  } else {
311  get_data_to_node = [this, th, &vals](EntityHandle tet) {
313  CHKERR mwlsMoab.tag_get_data(th, &tet, 1, &*vals.begin());
315  };
316  }
317 
318  // clean tag on nodes
319  CHKERR mwlsMoab.tag_clear_data(th, levelNodes, &*vals.begin());
320  for (auto tet : levelTets) {
321 
322  CHKERR get_data_to_node(tet);
323  int num_nodes;
324  const EntityHandle *conn;
325  CHKERR mwlsMoab.get_connectivity(tet, conn, num_nodes, false);
326  EntityHandle ho_node;
327  EntityType tet_type = mwlsMoab.type_from_handle(tet);
328  CHKERR mwlsMoab.high_order_node(tet, conn, tet_type, ho_node);
329  CHKERR mwlsMoab.tag_set_data(th, &ho_node, 1, &*vals.begin());
330  CHKERR mwlsMoab.tag_set_data(thMidNode, &tet, 1, &ho_node);
331  }
332 
333  } else {
334 
335  PetscBool add_analytical_internal_stress_operators = PETSC_FALSE;
336  CHKERR PetscOptionsGetBool(PETSC_NULL, "",
337  "-add_analytical_internal_stress_operators",
338  &add_analytical_internal_stress_operators, NULL);
339 
340  CHKERR mwlsMoab.tag_clear_data(th, levelNodes, &*vals.begin());
341  for (auto tet : levelTets) {
342  int num_nodes;
343  const EntityHandle *conn;
344  CHKERR mwlsMoab.get_connectivity(tet, conn, num_nodes, false);
345  EntityHandle ho_node;
346  EntityType tet_type = mwlsMoab.type_from_handle(tet);
347  CHKERR mwlsMoab.high_order_node(tet, conn, tet_type, ho_node);
348  CHKERR mwlsMoab.tag_set_data(thMidNode, &tet, 1, &ho_node);
349 
350  if (add_analytical_internal_stress_operators == PETSC_FALSE) {
351  VectorDouble coords(3);
352  CHKERR mwlsMoab.get_coords(&ho_node, 1, &coords[0]);
354  &coords[0], &coords[1], &coords[2]);
355  auto t_strain = CrackPropagation::analyticalStrainFunction(t_coords);
356 
357  constexpr double K = 166.666666666e9; // K = E / 3 / (1 - 2*nu),
358  // E = 200 GPa, nu = 0.3
359 
360  // Only diagonal, vectorial notation
361  vals[0] = 3 * K * t_strain(0, 0);
362  vals[1] = 3 * K * t_strain(1, 1);
363  vals[2] = 3 * K * t_strain(2, 2);
364  }
365 
366  CHKERR mwlsMoab.tag_set_data(th, &ho_node, 1, &*vals.begin());
367  }
368  }
369 
371 }

◆ loadMWLSMesh()

MoFEMErrorCode FractureMechanics::MWLSApprox::loadMWLSMesh ( std::string  file_name)

Load mesh with data and do basic calculation

Parameters
file_nameFile name
Returns
Error code

Definition at line 139 of file MWLS.cpp.

139  {
141 
143 
144  const char *option;
145  option = "";
146  CHKERR mwlsMoab.load_file(file_name.c_str(), 0, option);
147  CHKERR mwlsMoab.get_entities_by_dimension(0, 3, levelTets);
148  if (levelTets.empty()) {
150  "No 3d element in MWLS mesh");
151  }
152 
153  CHKERR mwlsMoab.get_adjacencies(levelTets, 1, true, levelEdges,
154  moab::Interface::UNION);
155 
156  // Create HO node on tetrahedral
157  if (!useNodalData) {
158  EntityHandle meshset;
159  CHKERR mwlsMoab.create_meshset(MESHSET_SET, meshset);
160  CHKERR mwlsMoab.add_entities(meshset, levelTets);
161  CHKERR mwlsMoab.convert_entities(meshset, false, false, true);
162  CHKERR mwlsMoab.delete_entities(&meshset, 1);
163 
164  // Get HO nodes in TET center
165  std::vector<double> edge_length;
166  for (auto tet : levelTets) {
167  int num_nodes;
168  const EntityHandle *conn;
169  CHKERR mwlsMoab.get_connectivity(tet, conn, num_nodes, false);
170  EntityHandle ho_node;
171  EntityType tit_type = mwlsMoab.type_from_handle(tet);
172  CHKERR mwlsMoab.high_order_node(tet, conn, tit_type, ho_node);
173  levelNodes.insert(ho_node);
174  }
175 
176  } else {
177  Range tets;
178  CHKERR mwlsMoab.get_entities_by_dimension(0, 3, tets);
179  CHKERR mwlsMoab.get_connectivity(tets, levelNodes, true);
180  }
181 
182  // Store edge nodes coordinates in FTensor
183  double edge_node_coords[6];
184  FTensor::Tensor1<double *, 3> t_node_edge[2] = {
185  FTensor::Tensor1<double *, 3>(edge_node_coords, &edge_node_coords[1],
186  &edge_node_coords[2]),
187  FTensor::Tensor1<double *, 3>(&edge_node_coords[3], &edge_node_coords[4],
188  &edge_node_coords[5])};
189 
191 
192  // Get edge lengths
193  maxEdgeL = 0;
194  for (auto edge : levelEdges) {
195  int num_nodes;
196  const EntityHandle *conn;
197  CHKERR mwlsMoab.get_connectivity(edge, conn, num_nodes, true);
198  CHKERR mwlsMoab.get_coords(conn, num_nodes, edge_node_coords);
199  t_node_edge[0](i) -= t_node_edge[1](i);
200  double l = sqrt(t_node_edge[0](i) * t_node_edge[0](i));
201  maxEdgeL = (maxEdgeL < l) ? l : maxEdgeL;
202  }
203 
204  // Create tree and find maximal edge length. Maximal edge length is used
205  // to find neighbours nodes to material point.
206  myTree = boost::make_shared<BVHTree>(&mwlsMoab);
207  {
208  std::ostringstream opts;
209  opts << "MAX_DEPTH=" << maxThreeDepth
210  // << ";MAX_PER_LEAF=" << maxElemsInDMFactor * nbBasePolynomials + 1
211  << ";SPLITS_PER_DIR=" << splitsPerDir;
212  FileOptions tree_opts(opts.str().c_str());
213 
214  MOFEM_LOG_CHANNEL("MWLSWorld");
215  MOFEM_LOG_TAG("MWLSWorld", "MWLS");
216  BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
217  MOFEM_LOG("MWLSWorld", Sev::verbose) << "Build tree ... ";
218  CHKERR myTree->build_tree(levelTets, &treeRoot, &tree_opts);
219  MOFEM_LOG("MWLSWorld", Sev::verbose) << "done";
220  }
221 
222  // moab::BoundBox box;
223  // CHKERR myTree->get_bounding_box(box);
224  // cout << box << endl;
225  // myTree->print();
226 
227  auto map_elems_positions = [&]() {
228  std::map<EntityHandle, std::array<double, 3>> elem_coords_map;
229  Range tets;
230  CHKERR mwlsMoab.get_entities_by_dimension(0, 3, tets, false);
231  MatrixDouble elem_coords;
232  for (auto p = tets.pair_begin(); p != tets.pair_end(); ++p) {
233  EntityHandle *conn;
234  int count, verts_per_e;
235  Range slot(p->first, p->second);
236  CHKERR mwlsMoab.connect_iterate(slot.begin(), slot.end(), conn,
237  verts_per_e, count);
238  if (count != (int)slot.size())
239  THROW_MESSAGE("Should be one continuos sequence");
240  elem_coords.resize(verts_per_e, 3, false);
241 
242  for (auto t = 0; t != count; ++t) {
243  CHKERR mwlsMoab.get_coords(conn, verts_per_e,
244  &*elem_coords.data().begin());
245  auto get_center = [&]() {
246  std::array<double, 3> c{0, 0, 0};
247  for (auto d : {0, 1, 2})
248  for (auto n = 0; n < verts_per_e; ++n)
249  c[d] += elem_coords(n, d);
250  for (auto d : {0, 1, 2})
251  c[d] /= verts_per_e;
252  return c;
253  };
254  elem_coords_map.emplace(std::make_pair(p->first + t, get_center()));
255  conn += verts_per_e;
256  }
257  }
258  return elem_coords_map;
259  };
260 
261  elemCentreMap = map_elems_positions();
262 
263  // Set tag for MWLS stress when SSLV116 test is run
265  Tag th;
266  std::array<double, 9> def;
267  def.fill(0);
268  CHKERR mwlsMoab.tag_get_handle("MED_SSLV116", 9, MB_TYPE_DOUBLE, th,
269  MB_TAG_CREAT | MB_TAG_SPARSE, def.data());
270  }
271 
273 }

Member Data Documentation

◆ A

ublas::symmetric_matrix<double, ublas::lower> FractureMechanics::MWLSApprox::A
private

Definition at line 641 of file MWLS.hpp.

◆ approxBasePoint

MatrixDouble FractureMechanics::MWLSApprox::approxBasePoint

Definition at line 213 of file MWLS.hpp.

◆ approxFun

VectorDouble FractureMechanics::MWLSApprox::approxFun
private

Definition at line 647 of file MWLS.hpp.

◆ approxPointCoords

VectorDouble3 FractureMechanics::MWLSApprox::approxPointCoords

Definition at line 214 of file MWLS.hpp.

◆ arcLengthDof

boost::weak_ptr<DofEntity> FractureMechanics::MWLSApprox::arcLengthDof

Definition at line 60 of file MWLS.hpp.

◆ B

MatrixDouble FractureMechanics::MWLSApprox::B
private

Definition at line 645 of file MWLS.hpp.

◆ baseFun

VectorDouble FractureMechanics::MWLSApprox::baseFun
private

Definition at line 646 of file MWLS.hpp.

◆ baseFunInvA

VectorDouble FractureMechanics::MWLSApprox::baseFunInvA
private

Definition at line 649 of file MWLS.hpp.

◆ diffA

ublas::symmetric_matrix<double, ublas::lower, ublas::row_major, DoubleAllocator> FractureMechanics::MWLSApprox::diffA[3]
private

Definition at line 653 of file MWLS.hpp.

◆ diffApproxFun

MatrixDouble FractureMechanics::MWLSApprox::diffApproxFun
private

Definition at line 660 of file MWLS.hpp.

◆ diffB

MatrixDouble FractureMechanics::MWLSApprox::diffB[3]
private

Definition at line 657 of file MWLS.hpp.

◆ diffBaseFun

VectorDouble FractureMechanics::MWLSApprox::diffBaseFun[3]
private

Definition at line 658 of file MWLS.hpp.

◆ diffDiffApproxFun

MatrixDouble FractureMechanics::MWLSApprox::diffDiffApproxFun
private

Definition at line 661 of file MWLS.hpp.

◆ diffDiffBaseFun

VectorDouble FractureMechanics::MWLSApprox::diffDiffBaseFun[9]
private

Definition at line 659 of file MWLS.hpp.

◆ diffDiffRhoAtGaussPts

boost::shared_ptr<MatrixDouble> FractureMechanics::MWLSApprox::diffDiffRhoAtGaussPts

Definition at line 208 of file MWLS.hpp.

◆ diffInvA

ublas::symmetric_matrix<double, ublas::lower, ublas::row_major, DoubleAllocator> FractureMechanics::MWLSApprox::diffInvA[3]
private

Definition at line 656 of file MWLS.hpp.

◆ diffRhoAtGaussPts

boost::shared_ptr<MatrixDouble> FractureMechanics::MWLSApprox::diffRhoAtGaussPts

Definition at line 207 of file MWLS.hpp.

◆ diffStressAtGaussPts

MatrixDouble FractureMechanics::MWLSApprox::diffStressAtGaussPts

Definition at line 203 of file MWLS.hpp.

◆ distLeafs

std::vector<std::pair<double, EntityHandle> > FractureMechanics::MWLSApprox::distLeafs
private

Definition at line 826 of file MWLS.hpp.

◆ dmFactor

double FractureMechanics::MWLSApprox::dmFactor
private

Relative value of weight function radius.

Definition at line 617 of file MWLS.hpp.

◆ dmNodesMap

std::map<EntityHandle, std::vector<double> > FractureMechanics::MWLSApprox::dmNodesMap

Store weight function radius at the nodes.

Definition at line 236 of file MWLS.hpp.

◆ elemCentreMap

std::map<EntityHandle, std::array<double, 3> > FractureMechanics::MWLSApprox::elemCentreMap
private

Definition at line 637 of file MWLS.hpp.

◆ F_lambda

Vec FractureMechanics::MWLSApprox::F_lambda

Definition at line 58 of file MWLS.hpp.

◆ functionTestHack

boost::function<VectorDouble(const double, const double, const double)> FractureMechanics::MWLSApprox::functionTestHack

This is used to set up analytical function for testing approximation.

Definition at line 602 of file MWLS.hpp.

◆ influenceNodes

std::vector<EntityHandle> FractureMechanics::MWLSApprox::influenceNodes
private

Definition at line 636 of file MWLS.hpp.

◆ influenceNodesData

MatrixDouble FractureMechanics::MWLSApprox::influenceNodesData
private

Definition at line 663 of file MWLS.hpp.

◆ influenceNodesMap

std::map<EntityHandle, std::vector<std::vector<EntityHandle> > > FractureMechanics::MWLSApprox::influenceNodesMap

Influence nodes on elements on at integration points.

Definition at line 230 of file MWLS.hpp.

◆ invA

MatrixDouble FractureMechanics::MWLSApprox::invA
private

Definition at line 643 of file MWLS.hpp.

◆ invAB

MatrixDouble FractureMechanics::MWLSApprox::invAB
private

Definition at line 648 of file MWLS.hpp.

◆ invABMap

std::map<EntityHandle, std::vector<MatrixDouble> > FractureMechanics::MWLSApprox::invABMap

Store Coefficient of base functions at integration points.

Definition at line 223 of file MWLS.hpp.

◆ L

ublas::triangular_matrix<double, ublas::lower> FractureMechanics::MWLSApprox::L
private

Definition at line 644 of file MWLS.hpp.

◆ leafsDist

std::vector<double> FractureMechanics::MWLSApprox::leafsDist
private

Definition at line 825 of file MWLS.hpp.

◆ leafsTetsCentre

std::vector<double> FractureMechanics::MWLSApprox::leafsTetsCentre
private

Definition at line 823 of file MWLS.hpp.

◆ leafsTetsVecLeaf

std::vector<EntityHandle> FractureMechanics::MWLSApprox::leafsTetsVecLeaf
private

Definition at line 822 of file MWLS.hpp.

◆ leafsVec

std::vector<EntityHandle> FractureMechanics::MWLSApprox::leafsVec
private

Definition at line 824 of file MWLS.hpp.

◆ levelEdges

Range FractureMechanics::MWLSApprox::levelEdges
private

Definition at line 635 of file MWLS.hpp.

◆ levelNodes

Range FractureMechanics::MWLSApprox::levelNodes
private

Definition at line 635 of file MWLS.hpp.

◆ levelTets

Range FractureMechanics::MWLSApprox::levelTets
private

Definition at line 635 of file MWLS.hpp.

◆ maxEdgeL

double FractureMechanics::MWLSApprox::maxEdgeL
private

Definition at line 633 of file MWLS.hpp.

◆ maxElemsInDMFactor

int FractureMechanics::MWLSApprox::maxElemsInDMFactor
private

Maximal number of elements in the horizon is calculated from max_nb_elems = maxElemsInDMFactor *nbBasePolynomials. Value is set by command line -mwls_max_elems_factor

Definition at line 620 of file MWLS.hpp.

◆ maxThreeDepth

int FractureMechanics::MWLSApprox::maxThreeDepth
private

Maximal three depths.

Definition at line 625 of file MWLS.hpp.

◆ mField

MoFEM::Interface& FractureMechanics::MWLSApprox::mField

Definition at line 55 of file MWLS.hpp.

◆ mwlsAnalyticalInternalStressTest

PetscBool FractureMechanics::MWLSApprox::mwlsAnalyticalInternalStressTest
private

Definition at line 630 of file MWLS.hpp.

◆ mwlsCore

moab::Core FractureMechanics::MWLSApprox::mwlsCore

Definition at line 56 of file MWLS.hpp.

◆ mwlsMaterialPositions

MatrixDouble FractureMechanics::MWLSApprox::mwlsMaterialPositions

Definition at line 212 of file MWLS.hpp.

◆ mwlsMoab

moab::Interface& FractureMechanics::MWLSApprox::mwlsMoab

Definition at line 57 of file MWLS.hpp.

◆ myTree

boost::shared_ptr<BVHTree> FractureMechanics::MWLSApprox::myTree
private

Definition at line 611 of file MWLS.hpp.

◆ nbBasePolynomials

int FractureMechanics::MWLSApprox::nbBasePolynomials
private

Number of base functions (1 - linear, 4 - linear, 10 - quadratic)

Definition at line 614 of file MWLS.hpp.

◆ nearestInfluenceNode

EntityHandle FractureMechanics::MWLSApprox::nearestInfluenceNode
private

Definition at line 639 of file MWLS.hpp.

◆ outData

VecVals FractureMechanics::MWLSApprox::outData
private

Definition at line 664 of file MWLS.hpp.

◆ outDeviatoricDifference

VecVals FractureMechanics::MWLSApprox::outDeviatoricDifference
private

Definition at line 668 of file MWLS.hpp.

◆ outDeviatoricError

double FractureMechanics::MWLSApprox::outDeviatoricError
private

Definition at line 669 of file MWLS.hpp.

◆ outDiffData

MatDiffVals FractureMechanics::MWLSApprox::outDiffData
private

Definition at line 671 of file MWLS.hpp.

◆ outDiffDiffData

MatDiffDiffVals FractureMechanics::MWLSApprox::outDiffDiffData
private

Definition at line 672 of file MWLS.hpp.

◆ outHydroStaticError

double FractureMechanics::MWLSApprox::outHydroStaticError
private

Definition at line 666 of file MWLS.hpp.

◆ rhoAtGaussPts

boost::shared_ptr<VectorDouble> FractureMechanics::MWLSApprox::rhoAtGaussPts

Definition at line 206 of file MWLS.hpp.

◆ shortenDm

double FractureMechanics::MWLSApprox::shortenDm
private

Radius of weight function used to calculate base function.

Definition at line 618 of file MWLS.hpp.

◆ singularCurrentDisplacement

MatrixDouble FractureMechanics::MWLSApprox::singularCurrentDisplacement

Definition at line 211 of file MWLS.hpp.

◆ singularInitialDisplacement

MatrixDouble FractureMechanics::MWLSApprox::singularInitialDisplacement

Definition at line 210 of file MWLS.hpp.

◆ splitsPerDir

int FractureMechanics::MWLSApprox::splitsPerDir
private

Split per direction in the tree.

Definition at line 626 of file MWLS.hpp.

◆ stressAtGaussPts

MatrixDouble FractureMechanics::MWLSApprox::stressAtGaussPts

Definition at line 202 of file MWLS.hpp.

◆ testA

ublas::symmetric_matrix<double, ublas::lower> FractureMechanics::MWLSApprox::testA
private

Definition at line 642 of file MWLS.hpp.

◆ tetsInBlock

Range FractureMechanics::MWLSApprox::tetsInBlock

Definition at line 215 of file MWLS.hpp.

◆ thMidNode

Tag FractureMechanics::MWLSApprox::thMidNode

Definition at line 217 of file MWLS.hpp.

◆ treeRoot

EntityHandle FractureMechanics::MWLSApprox::treeRoot
private

Definition at line 632 of file MWLS.hpp.

◆ treeTets

std::vector<EntityHandle> FractureMechanics::MWLSApprox::treeTets
private

Definition at line 821 of file MWLS.hpp.

◆ trisInBlock

Range FractureMechanics::MWLSApprox::trisInBlock

Definition at line 216 of file MWLS.hpp.

◆ useGlobalBaseAtMaterialReferenceConfiguration

PetscBool FractureMechanics::MWLSApprox::useGlobalBaseAtMaterialReferenceConfiguration
private

Definition at line 628 of file MWLS.hpp.

◆ useNodalData

PetscBool FractureMechanics::MWLSApprox::useNodalData
private

Definition at line 629 of file MWLS.hpp.


The documentation for this struct was generated from the following files:
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::K
VectorDouble K
Definition: Projection10NodeCoordsOnField.cpp:125
FractureMechanics::MWLSApprox::evalWeight
double evalWeight(const double r)
Definition: MWLS.hpp:781
FractureMechanics::MWLSApprox::evalDiffWeight
double evalDiffWeight(const double r)
Definition: MWLS.hpp:797
FractureMechanics::MWLSApprox::thMidNode
Tag thMidNode
Definition: MWLS.hpp:217
FractureMechanics::MWLSApprox::elemCentreMap
std::map< EntityHandle, std::array< double, 3 > > elemCentreMap
Definition: MWLS.hpp:637
FractureMechanics::MWLSApprox::levelNodes
Range levelNodes
Definition: MWLS.hpp:635
FractureMechanics::MWLSApprox::invAB
MatrixDouble invAB
Definition: MWLS.hpp:648
FTensor::Tensor1< double, 3 >
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
EntityHandle
FractureMechanics::MWLSApprox::nbBasePolynomials
int nbBasePolynomials
Definition: MWLS.hpp:614
young_modulus
double young_modulus
Young modulus.
Definition: plastic.cpp:121
FractureMechanics::MWLSApprox::diffB
MatrixDouble diffB[3]
Definition: MWLS.hpp:657
FractureMechanics::MWLSApprox::diffDiffApproxFun
MatrixDouble diffDiffApproxFun
Definition: MWLS.hpp:661
FractureMechanics::MWLSApprox::getMWLSOptions
MoFEMErrorCode getMWLSOptions()
get options from command line for MWLS
Definition: MWLS.cpp:71
FractureMechanics::MWLSApprox::F_lambda
Vec F_lambda
Definition: MWLS.hpp:58
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
FractureMechanics::MWLSApprox::B
MatrixDouble B
Definition: MWLS.hpp:645
FractureMechanics::MWLSApprox::myTree
boost::shared_ptr< BVHTree > myTree
Definition: MWLS.hpp:611
FractureMechanics::MWLSApprox::splitsPerDir
int splitsPerDir
Split per direction in the tree.
Definition: MWLS.hpp:626
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
FractureMechanics::MWLSApprox::diffRhoAtGaussPts
boost::shared_ptr< MatrixDouble > diffRhoAtGaussPts
Definition: MWLS.hpp:207
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
FractureMechanics::MWLSApprox::mwlsCore
moab::Core mwlsCore
Definition: MWLS.hpp:56
MoFEM::th
Tag th
Definition: Projection10NodeCoordsOnField.cpp:122
FractureMechanics::MWLSApprox::outData
VecVals outData
Definition: MWLS.hpp:664
FractureMechanics::MWLSApprox::calculateBase
MoFEMErrorCode calculateBase(const T &t_coords)
Definition: MWLS.hpp:675
FTensor::Tensor2_symmetric
Definition: Tensor2_symmetric_value.hpp:13
FractureMechanics::MWLSApprox::diffApproxFun
MatrixDouble diffApproxFun
Definition: MWLS.hpp:660
THROW_MESSAGE
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:574
FractureMechanics::MWLSApprox::leafsVec
std::vector< EntityHandle > leafsVec
Definition: MWLS.hpp:824
sdf.r
int r
Definition: sdf.py:8
FractureMechanics::MWLSApprox::leafsTetsVecLeaf
std::vector< EntityHandle > leafsTetsVecLeaf
Definition: MWLS.hpp:822
FractureMechanics::MWLSApprox::approxPointCoords
VectorDouble3 approxPointCoords
Definition: MWLS.hpp:214
FractureMechanics::MWLSApprox::outDiffData
MatDiffVals outDiffData
Definition: MWLS.hpp:671
FractureMechanics::MWLSApprox::L
ublas::triangular_matrix< double, ublas::lower > L
Definition: MWLS.hpp:644
FractureMechanics::MWLSApprox::influenceNodesData
MatrixDouble influenceNodesData
Definition: MWLS.hpp:663
FractureMechanics::MWLSApprox::influenceNodes
std::vector< EntityHandle > influenceNodes
Definition: MWLS.hpp:636
FractureMechanics::MWLSApprox::approxFun
VectorDouble approxFun
Definition: MWLS.hpp:647
FractureMechanics::MWLSApprox::levelEdges
Range levelEdges
Definition: MWLS.hpp:635
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
FractureMechanics::MWLSApprox::diffInvA
ublas::symmetric_matrix< double, ublas::lower, ublas::row_major, DoubleAllocator > diffInvA[3]
Definition: MWLS.hpp:656
FractureMechanics::MWLSApprox::rhoAtGaussPts
boost::shared_ptr< VectorDouble > rhoAtGaussPts
Definition: MWLS.hpp:206
a
constexpr double a
Definition: approx_sphere.cpp:30
FractureMechanics::MWLSApprox::arcLengthDof
boost::weak_ptr< DofEntity > arcLengthDof
Definition: MWLS.hpp:60
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
FractureMechanics::MWLSApprox::A
ublas::symmetric_matrix< double, ublas::lower > A
Definition: MWLS.hpp:641
FractureMechanics::MWLSApprox::leafsDist
std::vector< double > leafsDist
Definition: MWLS.hpp:825
FractureMechanics::MWLSApprox::dmFactor
double dmFactor
Relative value of weight function radius.
Definition: MWLS.hpp:617
poisson_ratio
double poisson_ratio
Poisson ratio.
Definition: plastic.cpp:122
FractureMechanics::MWLSApprox::maxThreeDepth
int maxThreeDepth
Maximal three depths.
Definition: MWLS.hpp:625
FractureMechanics::MWLSApprox::treeRoot
EntityHandle treeRoot
Definition: MWLS.hpp:632
FractureMechanics::MWLSApprox::mwlsMoab
moab::Interface & mwlsMoab
Definition: MWLS.hpp:57
FractureMechanics::MWLSApprox::testA
ublas::symmetric_matrix< double, ublas::lower > testA
Definition: MWLS.hpp:642
MoFEM::solve
static PetscErrorCode solve(Mat mat, Vec x, Vec y)
Definition: Schur.cpp:1223
FractureMechanics::MWLSApprox::invA
MatrixDouble invA
Definition: MWLS.hpp:643
MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
FractureMechanics::MWLSApprox::calculateDiffDiffBase
MoFEMErrorCode calculateDiffDiffBase(const T &t_coords, const double dm)
Definition: MWLS.hpp:750
FractureMechanics::MWLSApprox::outDiffDiffData
MatDiffDiffVals outDiffDiffData
Definition: MWLS.hpp:672
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
t
constexpr double t
plate stiffness
Definition: plate.cpp:58
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
t_kd
constexpr auto t_kd
Definition: free_surface.cpp:137
FractureMechanics::MWLSApprox::maxEdgeL
double maxEdgeL
Definition: MWLS.hpp:633
FractureMechanics::CrackPropagation::analyticalStrainFunction
static FTensor::Tensor2_symmetric< double, 3 > analyticalStrainFunction(FTensor::Tensor1< FTensor::PackPtr< double *, 1 >, 3 > &t_coords)
Definition: CrackPropagation.cpp:10030
FractureMechanics::MWLSApprox::baseFunInvA
VectorDouble baseFunInvA
Definition: MWLS.hpp:649
FTensor::Index< 'i', 3 >
FractureMechanics::MWLSApprox::baseFun
VectorDouble baseFun
Definition: MWLS.hpp:646
convert.n
n
Definition: convert.py:82
FractureMechanics::MWLSApprox::shortenDm
double shortenDm
Definition: MWLS.hpp:618
FractureMechanics::MWLSApprox::maxElemsInDMFactor
int maxElemsInDMFactor
Definition: MWLS.hpp:620
FractureMechanics::MWLSApprox::mwlsAnalyticalInternalStressTest
PetscBool mwlsAnalyticalInternalStressTest
Definition: MWLS.hpp:630
Range
FTensor::dd
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
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
FractureMechanics::MWLSApprox::mField
MoFEM::Interface & mField
Definition: MWLS.hpp:55
FractureMechanics::MWLSApprox::diffA
ublas::symmetric_matrix< double, ublas::lower, ublas::row_major, DoubleAllocator > diffA[3]
Definition: MWLS.hpp:653
FractureMechanics::MWLSApprox::leafsTetsCentre
std::vector< double > leafsTetsCentre
Definition: MWLS.hpp:823
FractureMechanics::MWLSApprox::diffBaseFun
VectorDouble diffBaseFun[3]
Definition: MWLS.hpp:658
FractureMechanics::MWLSApprox::useNodalData
PetscBool useNodalData
Definition: MWLS.hpp:629
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
FTensor::Ddg< double, 3, 3 >
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
FractureMechanics::MWLSApprox::calculateDiffBase
MoFEMErrorCode calculateDiffBase(const T &t_coords, const double dm)
Definition: MWLS.hpp:701
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
cholesky_decompose
size_t cholesky_decompose(const MATRIX &A, TRIA &L)
decompose the symmetric positive definit matrix A into product L L^T.
Definition: cholesky.hpp:52
cholesky_solve
void cholesky_solve(const TRIA &L, VEC &x, ublas::lower)
solve system L L^T x = b inplace
Definition: cholesky.hpp:221
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
sdf_wavy_2d.w
int w
Definition: sdf_wavy_2d.py:6
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition: Kronecker_Delta.hpp:49
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
FractureMechanics::MWLSApprox::useGlobalBaseAtMaterialReferenceConfiguration
PetscBool useGlobalBaseAtMaterialReferenceConfiguration
Definition: MWLS.hpp:628
FractureMechanics::MWLSApprox::levelTets
Range levelTets
Definition: MWLS.hpp:635
FractureMechanics::MWLSApprox::treeTets
std::vector< EntityHandle > treeTets
Definition: MWLS.hpp:821
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
FractureMechanics::MWLSApprox::nearestInfluenceNode
EntityHandle nearestInfluenceNode
Definition: MWLS.hpp:639
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
FractureMechanics::MWLSApprox::diffDiffRhoAtGaussPts
boost::shared_ptr< MatrixDouble > diffDiffRhoAtGaussPts
Definition: MWLS.hpp:208
FractureMechanics::MWLSApprox::diffDiffBaseFun
VectorDouble diffDiffBaseFun[9]
Definition: MWLS.hpp:659
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
FractureMechanics::MWLSApprox::distLeafs
std::vector< std::pair< double, EntityHandle > > distLeafs
Definition: MWLS.hpp:826
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182