v0.13.1
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< EntityHandle > influenceNodes
 
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< EntityHandle > treeTets
 
std::vector< EntityHandle > leafsTetsVecLeaf
 
std::vector< doubleleafsTetsCentre
 
std::vector< EntityHandle > leafsVec
 
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}
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:332
PetscBool useGlobalBaseAtMaterialReferenceConfiguration
Definition: MWLS.hpp:628
PetscBool mwlsAnalyticalInternalStressTest
Definition: MWLS.hpp:630
VectorDouble3 approxPointCoords
Definition: MWLS.hpp:214
boost::weak_ptr< DofEntity > arcLengthDof
Definition: MWLS.hpp:60
MoFEM::Interface & mField
Definition: MWLS.hpp:55
double dmFactor
Relative value of weight function radius.
Definition: MWLS.hpp:617
boost::shared_ptr< MatrixDouble > diffRhoAtGaussPts
Definition: MWLS.hpp:207
boost::shared_ptr< MatrixDouble > diffDiffRhoAtGaussPts
Definition: MWLS.hpp:208
moab::Interface & mwlsMoab
Definition: MWLS.hpp:57
boost::shared_ptr< VectorDouble > rhoAtGaussPts
Definition: MWLS.hpp:206

◆ ~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++) {
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 {
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}
void cholesky_solve(const TRIA &L, VEC &x, ublas::lower)
solve system L L^T x = b inplace
Definition: cholesky.hpp:221
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
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
double w(const double g, const double t)
const double r
rate factor
MatrixDouble diffB[3]
Definition: MWLS.hpp:657
std::vector< EntityHandle > influenceNodes
Definition: MWLS.hpp:636
ublas::symmetric_matrix< double, ublas::lower, ublas::row_major, DoubleAllocator > diffA[3]
Definition: MWLS.hpp:653
ublas::symmetric_matrix< double, ublas::lower > A
Definition: MWLS.hpp:641
double evalDiffWeight(const double r)
Definition: MWLS.hpp:797
double evalWeight(const double r)
Definition: MWLS.hpp:781
EntityHandle nearestInfluenceNode
Definition: MWLS.hpp:639
ublas::triangular_matrix< double, ublas::lower > L
Definition: MWLS.hpp:644
ublas::symmetric_matrix< double, ublas::lower > testA
Definition: MWLS.hpp:642
virtual MPI_Comm & get_comm() const =0

◆ calculateBase()

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

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 
)
private

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 }
VectorDouble diffBaseFun[3]
Definition: MWLS.hpp:658

◆ calculateDiffDiffBase()

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

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 }
VectorDouble diffDiffBaseFun[9]
Definition: MWLS.hpp:659

◆ evalDiffDiffWeight()

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

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)
private

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}
MoFEMErrorCode calculateBase(const T &t_coords)
Definition: MWLS.hpp:675

◆ 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}
constexpr double a
VectorDouble baseFunInvA
Definition: MWLS.hpp:649
ublas::symmetric_matrix< double, ublas::lower, ublas::row_major, DoubleAllocator > diffInvA[3]
Definition: MWLS.hpp:656
MatrixDouble diffApproxFun
Definition: MWLS.hpp:660
MoFEMErrorCode calculateDiffBase(const T &t_coords, const double dm)
Definition: MWLS.hpp:701

◆ 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}
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
FTensor::Index< 'n', SPACE_DIM > n
MatrixDouble diffDiffApproxFun
Definition: MWLS.hpp:661
MoFEMErrorCode calculateDiffDiffBase(const T &t_coords, const double dm)
Definition: MWLS.hpp:750

◆ evalWeight()

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

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}
MatrixDouble influenceNodesData
Definition: MWLS.hpp:663

◆ 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 {
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}
MatDiffDiffVals outDiffDiffData
Definition: MWLS.hpp:672

◆ 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
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}
Kronecker Delta class symmetric.
FTensor::Index< 'm', SPACE_DIM > m
constexpr auto t_kd
FTensor::Index< 'l', 3 > l
double young_modulus
Definition: plastic.cpp:96
constexpr double t
plate stiffness
Definition: plate.cpp:59

◆ getInfluenceNodes() [1/2]

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

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 )",
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}
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:304
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
std::vector< EntityHandle > leafsVec
Definition: MWLS.hpp:824
boost::shared_ptr< BVHTree > myTree
Definition: MWLS.hpp:611
std::vector< std::pair< double, EntityHandle > > distLeafs
Definition: MWLS.hpp:826
std::vector< double > leafsTetsCentre
Definition: MWLS.hpp:823
std::vector< EntityHandle > leafsTetsVecLeaf
Definition: MWLS.hpp:822
std::vector< double > leafsDist
Definition: MWLS.hpp:825
std::map< EntityHandle, std::array< double, 3 > > elemCentreMap
Definition: MWLS.hpp:637
std::vector< EntityHandle > treeTets
Definition: MWLS.hpp:821

◆ getInvAB()

MatrixDouble & FractureMechanics::MWLSApprox::getInvAB ( )

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", "",
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}
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
int maxThreeDepth
Maximal three depths.
Definition: MWLS.hpp:625
int splitsPerDir
Split per direction in the tree.
Definition: MWLS.hpp:626

◆ getShortenDm()

auto & FractureMechanics::MWLSApprox::getShortenDm ( )

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

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 {
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}
static Index< 'K', 3 > K
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
static FTensor::Tensor2_symmetric< double, 3 > analyticalStrainFunction(FTensor::Tensor1< FTensor::PackPtr< double *, 1 >, 3 > &t_coords)

◆ 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}
static Index< 'p', 3 > p
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:561
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:277
MoFEMErrorCode getMWLSOptions()
get options from command line for MWLS
Definition: MWLS.cpp:71

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: