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

Collaboration diagram for FractureMechanics::MWLSApprox:

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< FaceElementForcesAndSourcesCore >	OpMWLSCalculateBaseCoeffcientsForFaceAtGaussPts

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 VecVals &	getDataApprox ()

const MatDiffVals &	getDiffDataApprox ()

const MatDiffDiffVals &	getDiffDiffDataApprox ()

MatrixDouble &	getInvAB ()

std::vector< EntityHandle > &	getInfluenceNodes ()

auto &	getShortenDm ()

Public Attributes
MoFEM::Interface &	mField

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< double >	leafsTetsCentre

std::vector< EntityHandle >	leafsVec

std::vector< double >	leafsDist

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

typedef OpMWLSCalculateBaseCoeffcientsAtGaussPtsTmpl< VolumeElementForcesAndSourcesCore> FractureMechanics::MWLSApprox::OpMWLSCalculateBaseCoeffcientsAtGaussPts

Definition at line 296 of file MWLS.hpp.

◆ OpMWLSCalculateBaseCoeffcientsForFaceAtGaussPts

typedef OpMWLSCalculateBaseCoeffcientsAtGaussPtsTmpl< FaceElementForcesAndSourcesCore> FractureMechanics::MWLSApprox::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.

    : mField(m_field), mwlsMoab(mwlsCore), F_lambda(PETSC_NULL),
      arcLengthDof(arc_length_dof), nbBasePolynomials(1), dmFactor(1),
      maxElemsInDMFactor(2),
      useGlobalBaseAtMaterialReferenceConfiguration(PETSC_FALSE),
      useNodalData(PETSC_FALSE), mwlsAnalyticalInternalStressTest(PETSC_FALSE),
      approxPointCoords(3) {
 
  if (!LogManager::checkIfChannelExist("MWLSWorld")) {
    auto core_log = logging::core::get();
 
    core_log->add_sink(
        LogManager::createSink(LogManager::getStrmWorld(), "MWLSWorld"));
    core_log->add_sink(
        LogManager::createSink(LogManager::getStrmSync(), "MWLSSync"));
    core_log->add_sink(
        LogManager::createSink(LogManager::getStrmSelf(), "MWLSSelf"));
 
    LogManager::setLog("MWLSWorld");
    LogManager::setLog("MWLSSync");
    LogManager::setLog("MWLSSelf");
 
    MOFEM_LOG_TAG("MWLSWorld", "MWLS");
    MOFEM_LOG_TAG("MWLSSync", "MWLS");
    MOFEM_LOG_TAG("MWLSSelf", "MWLS");
  }
 
  MOFEM_LOG("MWLSWorld", Sev::noisy) << "MWLS created";
 
  rhoAtGaussPts = boost::make_shared<VectorDouble>();
  diffRhoAtGaussPts = boost::make_shared<MatrixDouble>();
  diffDiffRhoAtGaussPts = boost::make_shared<MatrixDouble>();
}

◆ ~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_coords material position

Returns: error code

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

Definition at line 576 of file MWLS.cpp.

                                                                         {
  MoFEMFunctionBegin;
 
  auto to_range = [&](auto ents) {
    Range r;
    r.insert_list(ents.begin(), ents.end());
    return r;
  };
 
  auto save_entities = [&](const std::string name, Range ents) {
    MoFEMFunctionBeginHot;
    if (!ents.empty()) {
      EntityHandle meshset;
      CHKERR mwlsMoab.create_meshset(MESHSET_SET, meshset);
      CHKERR mwlsMoab.add_entities(meshset, ents);
 
      CHKERR mwlsMoab.write_file(name.c_str(), "VTK", "", &meshset, 1);
      CHKERR mwlsMoab.delete_entities(&meshset, 1);
    }
    MoFEMFunctionReturnHot(0);
  };
 
  // Determine points in influence volume
  auto trim_nodes = [&](const double dm) {
    MoFEMFunctionBeginHot;
    FTensor::Index<'i', 3> i;
    FTensor::Tensor1<double, 3> t_approx_point(
        approxPointCoords[0], approxPointCoords[1], approxPointCoords[2]);
    const double dm2 = shortenDm * shortenDm;
    VectorDouble node_coors_vec(3 * influenceNodes.size());
    CHKERR mwlsMoab.get_coords(&*influenceNodes.begin(), influenceNodes.size(),
                               &*node_coors_vec.begin());
    FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_node_coords(
        &node_coors_vec[0], &node_coors_vec[1], &node_coors_vec[2]);
    decltype(influenceNodes) influence_nodes_tmp;
    influence_nodes_tmp.reserve(influenceNodes.size());
    for (auto node : influenceNodes) {
      t_node_coords(i) -= t_approx_point(i);
      const double r2 = t_node_coords(i) * t_node_coords(i);
      if (r2 < dm2)
        influence_nodes_tmp.emplace_back(node);
      ++t_node_coords;
    }
    influenceNodes.swap(influence_nodes_tmp);
    MoFEMFunctionReturnHot(0);
  };
 
  auto eval_AB = [this, global_derivatives](const double dm) {
    MoFEMFunctionBegin;
    FTensor::Index<'i', 3> i;
 
    A.resize(nbBasePolynomials, nbBasePolynomials, false);
    B.resize(nbBasePolynomials, influenceNodes.size(), false);
    A.clear();
    B.clear();
    for (int d = 0; d != 3; d++) {
      diffA[d].resize(nbBasePolynomials, nbBasePolynomials, false);
      diffB[d].resize(nbBasePolynomials, influenceNodes.size(), false);
      diffA[d].clear();
      diffB[d].clear();
    }
 
    FTensor::Tensor1<double, 3> t_approx_point(
        approxPointCoords[0], approxPointCoords[1], approxPointCoords[2]);
 
    VectorDouble nodes_coords(3 * influenceNodes.size());
    CHKERR mwlsMoab.get_coords(&*influenceNodes.begin(), influenceNodes.size(),
                               &*nodes_coords.begin());
    double *ptr = &*nodes_coords.begin();
    FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_node_coords(
        ptr, ptr + 1, ptr + 2);
 
    double min_distance = -1;
    nearestInfluenceNode = 0;
 
    for (int ii = 0; ii != influenceNodes.size(); ++ii) {
      t_node_coords(i) -= t_approx_point(i);
      t_node_coords(i) /= shortenDm;
      const double r = sqrt(t_node_coords(i) * t_node_coords(i));
 
      // Get the influence node that has minimum distance from the node in the
      // projected mesh
      if (!ii || r < min_distance) {
        min_distance = r;
        nearestInfluenceNode = influenceNodes[ii];
      }
 
      FTensor::Tensor1<double, 3> t_diff_r;
      if (r > 0) {
        t_diff_r(i) = (1. / r) * t_node_coords(i) * (-1 / shortenDm);
      } else {
        t_diff_r(i) = 0;
      }
      // Weights
      const double w = evalWeight(r);
      calculateBase<FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>>(
          t_node_coords);
 
      // iterate base functions
      // A_kj = sum_ii  w * b_k * b_j
      //
      // B_k = w * b_k
      //
      // diffA_kj = sum_ii diff_w * b_k * b_j + w * diff_b_k * b_j + w * b_k *
      // diff_b_j
      //
      // diffB =  diff_w * b_k + w * diff_bk
      for (int k = 0; k != nbBasePolynomials; k++) {
        const double wp = w * baseFun[k];
        for (int j = 0; j <= k; j++) {
          A(k, j) += wp * baseFun[j];
        }
        B(k, ii) = wp;
        if (global_derivatives) {
          const double diff_w_r = evalDiffWeight(r);
          const double diff_wp_r = diff_w_r * baseFun[k];
          for (int d = 0; d != 3; d++) {
            const double diff_wp_r_dx = diff_wp_r * t_diff_r(d);
            for (int j = 0; j <= k; j++) {
              diffA[d](k, j) += diff_wp_r_dx * baseFun[j];
            }
            diffB[d](k, ii) = diff_wp_r_dx;
          }
        }
      }
      ++t_node_coords;
    }
    if (nearestInfluenceNode == 0) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Nearest node not found \n");
    }
 
    MoFEMFunctionReturn(0);
  };
 
  // Invert matrix A
  auto invert_A = [this]() {
    MoFEMFunctionBegin;
 
    auto solve = [&](const size_t nb_base_functions) {
      A.resize(nb_base_functions, nb_base_functions, true);
      L.resize(nb_base_functions, nb_base_functions, false);
      if (cholesky_decompose(A, L)) {
        return false;
      } else {
 
        MatrixDouble inv_a(nb_base_functions, nb_base_functions, false);
        for (int k = 0; k != nb_base_functions; k++) {
          ublas::matrix_column<MatrixDouble> mr(inv_a, k);
          mr(k) = 1;
          cholesky_solve(L, mr, ublas::lower());
        }
 
        if (nb_base_functions == nbBasePolynomials) {
          invA.swap(inv_a);
        } else {
          invA.resize(nbBasePolynomials, nbBasePolynomials, false);
          invA.clear();
          for (size_t i = 0; i != nb_base_functions; ++i)
            for (size_t j = 0; j != nb_base_functions; ++j)
              invA(i, j) = inv_a(i, j);
        }
 
        return true;
      }
    };
 
    auto throw_error = [&]() {
      MoFEMFunctionBegin;
      FTensor::Tensor1<double, 3> t_approx_point(
          approxPointCoords[0], approxPointCoords[1], approxPointCoords[2]);
      cerr << "Point: " << t_approx_point << endl;
      cerr << "Matrix A: " << A << endl;
      SETERRQ1(mField.get_comm(), MOFEM_OPERATION_UNSUCCESSFUL,
               "Failed to invert matrix - solution could be "
               "to increase dmFactor to bigger value, e.g. -mwls_dm 4, "
               "or reduce of number of base fuctions "
               "-mwls_number_of_base_functions 1. The number of nodes "
               "found in "
               "radius is %u",
               influenceNodes.size());
      MoFEMFunctionReturn(0);
    };
 
    auto solve_one = [&]() {
      MoFEMFunctionBegin;
      if (!solve(1))
        CHKERR throw_error();
      MoFEMFunctionReturn(0);
    };
 
    testA = A;
    if (!solve(nbBasePolynomials)) {
      if (nbBasePolynomials == 10) {
        if (!solve(4))
          CHKERR solve_one();
      } else if (nbBasePolynomials == 4) {
        CHKERR solve_one();
      } else {
        CHKERR throw_error();
      }
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto calculate_shape_functions_coefficients = [this]() {
    MoFEMFunctionBegin;
    // Calculate base functions (shape functions) coefficients
    invAB.resize(invA.size1(), B.size2(), false);
    noalias(invAB) = prod(invA, B);
    MoFEMFunctionReturn(0);
  };
 
  // CHKERR trim_nodes(shortenDm);
  // CHKERR save_entities("tree_tets_trim.vtk", to_range(influenceNodes));
 
  CHKERR eval_AB(shortenDm);
  CHKERR invert_A();
  CHKERR calculate_shape_functions_coefficients();
 
  MoFEMFunctionReturn(0);
}

◆ 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.

                                                  {
    MoFEMFunctionBeginHot;
 
    baseFun.resize(nbBasePolynomials, false);
    baseFun.clear();
 
    baseFun[0] = 1;
    if (nbBasePolynomials > 1) {
      baseFun[1] = t_coords(0); // x
      baseFun[2] = t_coords(1); // y
      baseFun[3] = t_coords(2); // z
    }
 
    if (nbBasePolynomials > 4) {
      baseFun[4] = t_coords(0) * t_coords(0); // x^2
      baseFun[5] = t_coords(1) * t_coords(1); // y^2
      baseFun[6] = t_coords(2) * t_coords(2); // z^2
      baseFun[7] = t_coords(0) * t_coords(1); // xy
      baseFun[8] = t_coords(0) * t_coords(2); // xz
      baseFun[9] = t_coords(1) * t_coords(2); // yz
    }
 
    MoFEMFunctionReturnHot(0);
  }

◆ 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.

                                                                       {
    MoFEMFunctionBeginHot;
 
    // Local base functions
    for (int d = 0; d != 3; ++d) {
      diffBaseFun[d].resize(nbBasePolynomials, false);
      diffBaseFun[d].clear();
    }
 
    for (int d = 0; d != 3; d++) {
      diffBaseFun[d][0] = 0;
    }
 
    if (nbBasePolynomials > 1) {
      // derivative dx
      diffBaseFun[0][1] = 1. / dm;
      // derivative dy
      diffBaseFun[1][2] = 1. / dm;
      // derivative dz
      diffBaseFun[2][3] = 1. / dm;
    }
 
    if (nbBasePolynomials > 4) {
      // dx derivative
      diffBaseFun[0][4] = 2 * t_coords(0) / dm;
      diffBaseFun[0][5] = 0;
      diffBaseFun[0][6] = 0;
      diffBaseFun[0][7] = t_coords(1) / dm;
      diffBaseFun[0][8] = t_coords(2) / dm;
      diffBaseFun[0][9] = 0;
      // dy derivative
      diffBaseFun[1][4] = 0;
      diffBaseFun[1][5] = 2 * t_coords(1) / dm;
      diffBaseFun[1][6] = 0;
      diffBaseFun[1][7] = t_coords(0) / dm;
      diffBaseFun[1][8] = 0;
      diffBaseFun[1][9] = t_coords(2) / dm;
      // dz derivative
      diffBaseFun[2][4] = 0;
      diffBaseFun[2][5] = 0;
      diffBaseFun[2][6] = 2 * t_coords(2) / dm;
      diffBaseFun[2][7] = 0;
      diffBaseFun[2][8] = t_coords(0) / dm;
      diffBaseFun[2][9] = t_coords(1) / dm;
    }
    MoFEMFunctionReturnHot(0);
  }

◆ 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.

                                                                           {
    MoFEMFunctionBeginHot;
 
    // Local base functions
    for (int d = 0; d != 9; ++d) {
      diffDiffBaseFun[d].resize(nbBasePolynomials, false);
      diffDiffBaseFun[d].clear();
    }
 
    // if (nbBasePolynomials > 1) {
    // }
 
    if (nbBasePolynomials > 4) {
      const double dm2 = dm * dm;
      // set the non-zero values
      // dxd? derivative
      diffDiffBaseFun[0][4] = 2. / dm2; // dxdx
      diffDiffBaseFun[1][7] = 1. / dm2; // dxdy
      diffDiffBaseFun[2][8] = 1. / dm2; // dxdz
      // dyd? derivative
      diffDiffBaseFun[4][5] = 2. / dm2; // dydy
      diffDiffBaseFun[3][7] = 1. / dm2; // dydx
      diffDiffBaseFun[5][9] = 1. / dm2; // dydz
      // dzd? derivative
      diffDiffBaseFun[8][6] = 2. / dm2; // dzdz
      diffDiffBaseFun[6][8] = 1. / dm2; // dzdx
      diffDiffBaseFun[7][9] = 1. / dm2; // dzdy
    }
    MoFEMFunctionReturnHot(0);
  }

◆ evalDiffDiffWeight()

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

inlineprivate

Definition at line 812 of file MWLS.hpp.

                                                   {
    const double rr = r * r;
    if (r < 1) {
      return -36 * rr + 48 * r - 12;
    } else {
      return 0;
    }
  }

◆ evalDiffWeight()

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

inlineprivate

Definition at line 797 of file MWLS.hpp.

                                               {
    const double rr = r * r;
    const double rrr = rr * r;
    if (r < 1) {
      return -12 * r + 24 * rr - 12 * rrr;
    } else {
      return 0;
    }
    // if(r<=0.5) {
    //   return -8*r+12*rr;
    // } else if(r<=1) {
    //   return -4+8*r-4*rr;
    // } else return 0;
  }

◆ evaluateApproxFun()

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

evaluate approximation function at material point

Parameters

material_coords material position

Returns: error code

Definition at line 801 of file MWLS.cpp.

                                                                 {
  MoFEMFunctionBegin;
  // Determine points in influence volume
 
  auto cal_base_functions_at_eval_point = [this, eval_point](const double dm) {
    MoFEMFunctionBegin;
 
    FTensor::Tensor1<double, 3> t_approx_point(
        approxPointCoords[0], approxPointCoords[1], approxPointCoords[2]);
    FTensor::Tensor1<double, 3> t_eval_point(eval_point[0], eval_point[1],
                                             eval_point[2]);
 
    FTensor::Index<'i', 3> i;
    FTensor::Tensor1<double, 3> t_delta;
    t_delta(i) = t_eval_point(i) - t_approx_point(i);
    t_delta(i) /= dm;
 
    CHKERR calculateBase(t_delta);
    MoFEMFunctionReturn(0);
  };
 
  auto calculate_shape_functions = [this]() {
    MoFEMFunctionBegin;
    // Calculate approximation functions at eval point
    approxFun.resize(influenceNodes.size(), false);
    noalias(approxFun) = prod(baseFun, invAB);
    MoFEMFunctionReturn(0);
  };
 
  CHKERR cal_base_functions_at_eval_point(shortenDm);
  CHKERR calculate_shape_functions();
 
  MoFEMFunctionReturn(0);
}

◆ evaluateFirstDiffApproxFun()

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

evaluate 1st derivative of approximation function at material point

Parameters

global_derivatives decide whether global derivative will be calculated

Returns: error code

Definition at line 836 of file MWLS.cpp.

                                                                               {
  MoFEMFunctionBegin;
 
  auto cal_diff_base_functions_at_eval_point = [this,
                                                eval_point](const double dm) {
    MoFEMFunctionBegin;
 
    FTensor::Tensor1<double, 3> t_approx_point(
        approxPointCoords[0], approxPointCoords[1], approxPointCoords[2]);
    FTensor::Tensor1<double, 3> t_eval_point(eval_point[0], eval_point[1],
                                             eval_point[2]);
 
    FTensor::Index<'i', 3> i;
    FTensor::Tensor1<double, 3> t_delta;
    t_delta(i) = t_eval_point(i) - t_approx_point(i);
    t_delta(i) /= dm;
 
    CHKERR calculateDiffBase(t_delta, dm);
    MoFEMFunctionReturn(0);
  };
 
  auto calculate_global_shape_functions_derivatives = [this]() {
    MoFEMFunctionBegin;
 
    // Calculate derivatives of base functions
    VectorDouble tmp;
    tmp.resize(nbBasePolynomials, false);
    baseFunInvA.resize(invA.size2(), false);
    noalias(baseFunInvA) = prod(baseFun, invA);
 
    for (int d = 0; d != 3; d++) {
      diffInvA[d].resize(nbBasePolynomials, nbBasePolynomials, false);
      diffApproxFun.resize(3, influenceNodes.size(), false);
    }
 
    for (int d = 0; d != 3; d++) {
      MatrixDouble a = prod(diffA[d], invA);
      noalias(diffInvA[d]) = -prod(invA, a);
    }
 
    // Calculate line derivatives for each coordinate direction
    for (int d = 0; d != 3; ++d) {
 
      ublas::matrix_row<MatrixDouble> mr(diffApproxFun, d);
 
      // b,x * A^-1 B
      noalias(mr) = prod(diffBaseFun[d], invAB);
 
      // b * (A^-1 A,x A^-1) * B
      noalias(tmp) = prod(baseFun, diffInvA[d]);
      mr += prod(tmp, B);
 
      // b * A^-1 * B,x
      mr += prod(baseFunInvA, diffB[d]);
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto calculate_local_shape_functions_derivatives = [this]() {
    MoFEMFunctionBegin;
 
    // Calculate derivatives of base functions
    diffApproxFun.resize(3, influenceNodes.size(), false);
    // Calculate line derivatives for each coordinate direction
    for (int d = 0; d != 3; ++d) {
      ublas::matrix_row<MatrixDouble> mr(diffApproxFun, d);
      // b,x * A^-1 B
      noalias(mr) = prod(diffBaseFun[d], invAB);
    }
 
    MoFEMFunctionReturn(0);
  };
 
  CHKERR cal_diff_base_functions_at_eval_point(shortenDm);
  if (global_derivatives)
    CHKERR calculate_global_shape_functions_derivatives();
  else
    CHKERR calculate_local_shape_functions_derivatives();
 
  MoFEMFunctionReturn(0);
}

◆ evaluateSecondDiffApproxFun()

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

evaluate 2nd derivative of approximation function at material point

Parameters

global_derivatives decide whether global derivative will be calculated

Returns: error code

Definition at line 921 of file MWLS.cpp.

                                                                 {
  MoFEMFunctionBegin;
 
  auto cal_diff_diff_base_functions_at_eval_point =
      [this, eval_point](const double dm) {
        MoFEMFunctionBegin;
 
        FTensor::Tensor1<double, 3> t_approx_point(
            approxPointCoords[0], approxPointCoords[1], approxPointCoords[2]);
        FTensor::Tensor1<double, 3> t_eval_point(eval_point[0], eval_point[1],
                                                 eval_point[2]);
 
        FTensor::Index<'i', 3> i;
        FTensor::Tensor1<double, 3> t_delta;
        t_delta(i) = t_eval_point(i) - t_approx_point(i);
        t_delta(i) /= dm;
 
        CHKERR calculateDiffDiffBase(t_delta, dm);
        MoFEMFunctionReturn(0);
      };
 
  auto calculate_global_shape_functions_second_derivatives = [this]() {
    MoFEMFunctionBegin;
    // to be implemented
    MoFEMFunctionReturn(0);
  };
 
  auto calculate_local_shape_functions_second_derivatives = [this]() {
    MoFEMFunctionBegin;
 
    // Calculate derivatives of base functions
    // for (int d = 0; d != 3; d++) {
    diffDiffApproxFun.resize(9, influenceNodes.size(), false);
    // }
    // Calculate line derivatives for each coordinate direction
    for (int n = 0; n != 3; ++n) {
      for (int d = 0; d != 3; ++d) {
        const int indx = d + 3 * n;
        ublas::matrix_row<MatrixDouble> mr(diffDiffApproxFun, indx);
        // b,xy * A^-1 B
        noalias(mr) = prod(diffDiffBaseFun[indx], invAB);
      }
    }
    MoFEMFunctionReturn(0);
  };
 
  CHKERR cal_diff_diff_base_functions_at_eval_point(shortenDm);
  if (global_derivatives) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented yet");
    // CHKERR calculate_global_shape_functions_second_derivatives();
  } else
    CHKERR calculate_local_shape_functions_second_derivatives();
 
  MoFEMFunctionReturn(0);
}

◆ evalWeight()

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

inlineprivate

Definition at line 781 of file MWLS.hpp.

                                           {
    const double rr = r * r;
    const double rrr = rr * r;
    const double r4 = rrr * r;
    if (r < 1) {
      return 1 - 6 * rr + 8 * rrr - 3 * r4;
    } else {
      return 0;
    }
    // if(r<=0.5) {
    //   return 2./3-4*rr+4*rrr;
    // } else if(r<=1) {
    //   return 4./3-4*r+4*rr-(4./3)*rrr;
    // } else return 0;
  }

◆ getDataApprox()

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

Get values at point

Returns: error code

Definition at line 992 of file MWLS.cpp.

                                                   {
  noalias(outData) = prod(approxFun, influenceNodesData);
  return outData;
}

◆ getDiffDataApprox()

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

Get derivatives of values at points

Returns: error code

Definition at line 997 of file MWLS.cpp.

                                                           {
  noalias(outDiffData) = prod(diffApproxFun, influenceNodesData);
  return outDiffData;
}

◆ 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.

                                                                   {
  noalias(outDiffDiffData) = prod(diffDiffApproxFun, influenceNodesData);
  return outDiffDiffData;
}

◆ 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.

                                                            {
  MoFEMFunctionBegin;
 
  int length;
  CHKERR mwlsMoab.tag_get_length(th, length);
 
  FTensor::Index<'i', 3> i;
  FTensor::Index<'j', 3> j;
  FTensor::Index<'k', 3> k;
  FTensor::Index<'l', 3> l;
 
  auto get_compliance_matrix = [&]() {
    FTensor::Ddg<double, 3, 3> t_C;
    t_C(i, j, k, l) = 0;
 
    const double c = poisson_ratio / young_modulus;
    t_C(0, 0, 0, 0) = 1. / young_modulus;
    t_C(0, 0, 1, 1) = -c;
    t_C(0, 0, 2, 2) = -c;
 
    t_C(1, 1, 0, 0) = -c;
    t_C(1, 1, 1, 1) = 1. / young_modulus;
    t_C(1, 1, 2, 2) = -c;
 
    t_C(2, 2, 0, 0) = -c;
    t_C(2, 2, 1, 1) = -c;
    t_C(2, 2, 2, 2) = 1. / young_modulus;
 
    const double d = (1 + poisson_ratio) / young_modulus;
    t_C(0, 1, 0, 1) = d;
    t_C(0, 2, 0, 2) = d;
    t_C(1, 2, 1, 2) = d;
 
    return t_C;
  };
 
  auto t_C = get_compliance_matrix();
 
  if (length != outData.size())
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Tag length is not consistent with outData size. \nMost probably "
            "getTagData has not been invoked yet!\n");
 
  VectorDouble val_at_nearest_influence_node;
  val_at_nearest_influence_node.resize(length, false);
  val_at_nearest_influence_node.clear();
  CHKERR mwlsMoab.tag_get_data(th, &nearestInfluenceNode, 1,
                               &*val_at_nearest_influence_node.begin());
 
  // SIXX 0 (00) SIYY 1 (11) SIZZ 2 (22) SIXY 3 (01) SIXZ 4 (02) SIYZ 5 (12)
  // T* d00, T* d01, T* d02, T* d11, T* d12, T* d22
  auto get_symm_tensor = [](auto &m) {
    FTensor::Tensor2_symmetric<double *, 3> t(
 
        &m(0), &m(3), &m(4),
 
        &m(1), &m(5),
 
        &m(2)
 
    );
    return t;
  };
 
  auto t_total_stress_error = FTensor::Tensor2_symmetric<double, 3>();
  auto t_stress_at_nearset_node =
      get_symm_tensor(val_at_nearest_influence_node);
  auto t_mwls_stress = get_symm_tensor(outData);
 
  auto get_energy_norm = [&](auto &t_s) {
    return 0.5 * t_s(i, j) * (t_C(i, j, k, l) * t_s(k, l));
  };
  auto get_stress_nrm2 = [i, j](auto &t_s) { return t_s(i, j) * t_s(i, j); };
 
  total_stress_at_node = get_stress_nrm2(t_stress_at_nearset_node);
  total_energy = get_energy_norm(t_stress_at_nearset_node);
 
  t_total_stress_error(i, j) =
      t_mwls_stress(i, j) - t_stress_at_nearset_node(i, j);
  total_stress_error_at_node = get_stress_nrm2(t_total_stress_error);
  total_energy_error = get_energy_norm(t_total_stress_error);
 
  constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
  // calculate error of the hydro static part of the stress
  double hydro_static_difference =
      (t_kd(i, j) * t_total_stress_error(i, j)) / 3.;
  hydro_static_error_at_node =
      hydro_static_difference * hydro_static_difference;
 
  // calculate error of the deviator part of the stress
  auto t_dev_error = FTensor::Tensor2_symmetric<double, 3>();
  t_dev_error(i, j) =
      t_total_stress_error(i, j) - hydro_static_difference * t_kd(i, j);
  deviatoric_error_at_node = get_stress_nrm2(t_dev_error);
 
  MoFEMFunctionReturn(0);
}

◆ getInfluenceNodes() [1/2]

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

inline

Definition at line 605 of file MWLS.hpp.

                                                      {
    return influenceNodes;
  }

◆ getInfluenceNodes() [2/2]

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

Get node in influence domain.

Parameters

material_coords material position

Returns: error code

Definition at line 373 of file MWLS.cpp.

                                                                      {
  constexpr double eps_overlap = 0.01;
  const size_t max_nb_elems = maxElemsInDMFactor * nbBasePolynomials + 1;
  MoFEMFunctionBegin;
 
    using DistEntPair = std::pair<double, EntityHandle>;
    using DistEntMap = multi_index_container<
        DistEntPair,
        indexed_by<
 
            ordered_non_unique<member<DistEntPair, double, &DistEntPair::first>>
 
            >>;
    DistEntMap dist_map;
 
    // FIXME: That is potentisl source of problems, should be sparate
    // function for setting approximation base point, or simply should be set
    // outside of that function. Other functions will behave depending how
    // approxPointCoords is set. And here is somehow hidden from the view,
    // i.e. easy to overlook this. TODO: create a new function like:
    // MoFEMErrorCode setApproxBasePt(double material_coords[3]) ?
 
    for (int dd : {0, 1, 2})
      approxPointCoords[dd] = material_coords[dd];
 
    if (myTree) {
      shortenDm = maxEdgeL * dmFactor; // Influence radius
 
      auto find_leafs = [&](const auto dm) {
        leafsVec.clear();
        leafsDist.clear();
        CHKERR myTree->distance_search(material_coords, shortenDm, leafsVec,
                                       1.0e-10, 1.0e-6, &leafsDist);
        distLeafs.clear();
        distLeafs.reserve(leafsDist.size());
        for (size_t i = 0; i != leafsVec.size(); ++i) {
          distLeafs.emplace_back(leafsDist[i], leafsVec[i]);
        }
        std::sort(
            distLeafs.begin(), distLeafs.end(),
            [](const auto &a, const auto &b) { return a.first < b.first; });
        return distLeafs;
      };
 
      auto get_influence_nodes = [&](const auto &tets) {
        std::vector<EntityHandle> influence_nodes;
        if (!useNodalData) {
          influence_nodes.resize(tets.size());
          CHKERR mwlsMoab.tag_get_data(thMidNode, &*tets.begin(), tets.size(),
                                       &*influence_nodes.begin());
        } else {
          Range nodes;
          CHKERR mwlsMoab.get_connectivity(&*tets.begin(), tets.size(), nodes,
                                           true);
          influence_nodes.clear();
          influence_nodes.reserve(nodes.size());
          for (auto n : nodes)
            influence_nodes.emplace_back(n);
        }
        return influence_nodes;
      };
 
      auto to_range = [&](auto ents) {
        Range r;
        r.insert_list(ents.begin(), ents.end());
        return r;
      };
 
      auto save_entities = [&](const std::string name, Range ents) {
        MoFEMFunctionBeginHot;
        if (!ents.empty()) {
          EntityHandle meshset;
          CHKERR mwlsMoab.create_meshset(MESHSET_SET, meshset);
          CHKERR mwlsMoab.add_entities(meshset, ents);
 
          CHKERR mwlsMoab.write_file(name.c_str(), "VTK", "", &meshset, 1);
          CHKERR mwlsMoab.delete_entities(&meshset, 1);
        }
        MoFEMFunctionReturnHot(0);
      };
 
      auto leafs = find_leafs(shortenDm);
 
      treeTets.clear();
      std::vector<EntityHandle> leafs_tets;
      moab::BoundBox box;
 
      auto get_dm2_from_influence_points = [&]() {
        FTensor::Index<'i', 3> i;
        VectorDouble node_coors_vec(3 * influenceNodes.size());
        CHKERR mwlsMoab.get_coords(&*influenceNodes.begin(),
                                   influenceNodes.size(),
                                   &*node_coors_vec.begin());
        FTensor::Tensor1<double, 3> t_approx_point(
            material_coords[0], material_coords[1], material_coords[2]);
        FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_node_coords(
            &node_coors_vec[0], &node_coors_vec[1], &node_coors_vec[2]);
        double dm2 = 0;
        for (auto node : influenceNodes) {
          t_node_coords(i) -= t_approx_point(i);
          const double r2 = t_node_coords(i) * t_node_coords(i);
          if (dm2 < r2)
            dm2 = r2;
          ++t_node_coords;
        }
        return dm2;
      };
 
      const double shorten_dm2 = shortenDm * shortenDm;
      double set_dm2 = shorten_dm2;
      auto min_max_dist_it = dist_map.get<0>().begin();
      FTensor::Index<'i', 3> i;
 
      for (auto l : leafs) {
        if (dist_map.size() < max_nb_elems ||
            ((l.first * l.first) <
             (set_dm2 + std::numeric_limits<double>::epsilon()))) {
 
          leafsTetsVecLeaf.clear();
          CHKERR mwlsMoab.get_entities_by_dimension(l.second, 3,
                                                    leafsTetsVecLeaf, false);
 
          leafsTetsCentre.resize(3 * leafsTetsVecLeaf.size());
          {
            FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3> t_c{
                &leafsTetsCentre[0], &leafsTetsCentre[leafsTetsVecLeaf.size()],
                &leafsTetsCentre[2 * leafsTetsVecLeaf.size()]};
 
            for (auto tet : leafsTetsVecLeaf) {
              const auto &c = elemCentreMap[tet];
              for (auto ii : {0, 1, 2}) {
                t_c(ii) = c[ii] - material_coords[ii];
              }
              ++t_c;
            }
          }
 
          {
            FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3> t_c{
                &leafsTetsCentre[0], &leafsTetsCentre[leafsTetsVecLeaf.size()],
                &leafsTetsCentre[2 * leafsTetsVecLeaf.size()]};
 
            for (auto tet : leafsTetsVecLeaf) {
 
              if (dist_map.size() > max_nb_elems)
                set_dm2 = std::min(min_max_dist_it->first, shorten_dm2);
 
              if (t_c(0) < set_dm2 && t_c(1) < set_dm2 && t_c(2) < set_dm2) {
                double dm2 = t_c(i) * t_c(i);
                if (dm2 < set_dm2) {
                  min_max_dist_it =
                      dist_map.get<0>().emplace(std::make_pair(dm2, tet)).first;
                  if (dist_map.size() > max_nb_elems) {
                    int nb = std::distance(dist_map.get<0>().begin(),
                                           min_max_dist_it);
                    for (; nb < max_nb_elems; ++nb)
                      ++min_max_dist_it;
                  }
                }
              }
 
              ++t_c;
            }
          }
        } else {
          break;
        }
      }
 
      treeTets.clear();
      treeTets.reserve(dist_map.size());
      const double dm2 = set_dm2 + std::numeric_limits<double>::epsilon();
      auto hi_it = dist_map.get<0>().upper_bound(dm2);
      for (auto it = dist_map.get<0>().begin(); it != hi_it; ++it)
        treeTets.emplace_back(it->second);
 
      if (treeTets.empty()) {
        for (auto &it : dist_map)
          MOFEM_LOG("MWLSWorld", Sev::error)
              << "dm map " << it.first << " " << it.second;
 
        MOFEM_LOG("MWLSWorld", Sev::error) << "leafs found " << leafs.size();
        MOFEM_LOG("MWLSWorld", Sev::error)
            << "dist map size " << dist_map.size();
        MOFEM_LOG("MWLSWorld", Sev::error) << "dm2 " << dm2;
        MOFEM_LOG("MWLSWorld", Sev::error) << "shorten dm is " << shortenDm;
        MOFEM_LOG_C("MWLSWorld", Sev::error, "point:  ( %g %g %g )",
                    approxPointCoords(0), approxPointCoords(1),
                    approxPointCoords(2));
        SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No leafs found.");
      }
 
      influenceNodes = get_influence_nodes(treeTets);
      const auto dm2_from_influence_pts = get_dm2_from_influence_points();
      shortenDm = (1 + eps_overlap) * sqrt(dm2_from_influence_pts);
 
      // CHKERR save_entities("tree_tets.vtk", to_range(influenceNodes));
  } else {
    SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "kd-three not build");
  }
  MoFEMFunctionReturn(0);
}

◆ 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.

                                          {
  MoFEMFunctionBegin;
  ierr = PetscOptionsBegin(mField.get_comm(), "",
                           "Get options for moving least square approximation",
                           "none");
  CHKERRQ(ierr);
  CHKERR PetscOptionsScalar("-mwls_dm",
                            "edge length scaling for influence radius", "",
                            dmFactor, &dmFactor, PETSC_NULL);
  MOFEM_LOG_C("MWLSWorld", Sev::inform, "### Input parameter: -mwls_dm %6.4e",
              dmFactor);
 
  CHKERR PetscOptionsInt(
      "-mwls_number_of_base_functions",
      "1 = constant, 4 = linear, 10 = quadratic approximation", "",
      nbBasePolynomials, &nbBasePolynomials, PETSC_NULL);
  MOFEM_LOG_C("MWLSWorld", Sev::inform,
              "### Input parameter: -mwls_number_of_base_functions %d",
              nbBasePolynomials);
 
  CHKERR PetscOptionsInt("-mwls_max_elems_factor",
                         "Set max elements factor max_nb_elems = "
                         "maxElemsInDMFactor * nbBasePolynomials",
                         "", maxElemsInDMFactor, &maxElemsInDMFactor,
                         PETSC_NULL);
  MOFEM_LOG_C("MWLSWorld", Sev::inform,
              "### Input parameter: -mwls_max_elems_factor %d",
              maxElemsInDMFactor);
 
  CHKERR PetscOptionsBool(
      "-mwls_use_global_base_at_reference_configuration",
      "true local mwls base functions at reference configuration are used", "",
      useGlobalBaseAtMaterialReferenceConfiguration,
      &useGlobalBaseAtMaterialReferenceConfiguration, NULL);
  CHKERRQ(ierr);
 
  CHKERR PetscOptionsBool(
      "-mwls_use_nodal_data",
      "true nodal data values are used (without averaging to the midpoint)", "",
      useNodalData, &useNodalData, NULL);
  MOFEM_LOG_C("MWLSWorld", Sev::verbose,
              "### Input parameter: -mwls_use_nodal_data %d", useNodalData);
 
  CHKERR PetscOptionsBool(
      "-mwls_analytical_internal_stress",
      "use analytical strain field for internal stress mwls test", "",
      mwlsAnalyticalInternalStressTest, &mwlsAnalyticalInternalStressTest,
      NULL);
  MOFEM_LOG_C("MWLSWorld", Sev::verbose,
              "### Input parameter: -test_mwls_internal_stress_analytical %d",
              mwlsAnalyticalInternalStressTest);
 
  maxThreeDepth = 30;
  CHKERR PetscOptionsInt("-mwls_tree_depth", "maximal three depths", "",
                         maxThreeDepth, &maxThreeDepth, PETSC_NULL);
  MOFEM_LOG("MWLSWorld", Sev::verbose)
      << "Maximal MWLS three depths " << maxThreeDepth;
  splitsPerDir = 12;
  CHKERR PetscOptionsInt("-mwls_splits_per_dir", "splits per direction", "",
                         splitsPerDir, &splitsPerDir, PETSC_NULL);
  MOFEM_LOG("MWLSWorld", Sev::verbose)
      << "Splits per direction " << splitsPerDir;
 
  ierr = PetscOptionsEnd();
  CHKERRQ(ierr);
  MoFEMFunctionReturn(0);
}

◆ 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

th tag

Returns: error code

Definition at line 978 of file MWLS.cpp.

                                            {
  MoFEMFunctionBegin;
  int length;
  CHKERR mwlsMoab.tag_get_length(th, length);
  influenceNodesData.resize(influenceNodes.size(), length, false);
  CHKERR mwlsMoab.tag_get_data(th, &*influenceNodes.begin(),
                               influenceNodes.size(),
                               &influenceNodesData(0, 0));
  outData.resize(length, false);
  outDiffData.resize(3, length, false);
  outDiffDiffData.resize(9, length, false);
  MoFEMFunctionReturn(0);
}

◆ 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.

                                                                {
    return (useGlobalBaseAtMaterialReferenceConfiguration == PETSC_TRUE);
  }

◆ getValuesToNodes()

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

get values form elements to nodes

Parameters

th	tag with approximated data
th	tag with approximated data

Returns: error code

Definition at line 275 of file MWLS.cpp.

                                                  {
  MoFEMFunctionBegin;
  if (useNodalData)
    MoFEMFunctionReturnHot(0);
 
  CHKERR mwlsMoab.tag_get_handle("MID_NODE", 1, MB_TYPE_HANDLE, thMidNode,
                                 MB_TAG_CREAT | MB_TAG_DENSE);
 
  int length;
  CHKERR mwlsMoab.tag_get_length(th, length);
 
  VectorDouble vals(length);
  vals.clear();
 
  if (mwlsAnalyticalInternalStressTest == PETSC_FALSE) {
 
    string name;
    CHKERR mwlsMoab.tag_get_name(th, name);
 
    std::function<MoFEMErrorCode(EntityHandle)> get_data_to_node;
    if (name.compare("RHO") == 0) {
      get_data_to_node = [this, th, &vals](EntityHandle tet) {
        MoFEMFunctionBeginHot;
        Range nodes;
        double avg_val = 0;
        CHKERR mwlsMoab.get_connectivity(&tet, 1, nodes, true);
        for (auto &n : nodes) {
          double tag_data;
          CHKERR mwlsMoab.tag_get_data(th, &n, 1, &tag_data);
          avg_val += tag_data;
        }
        avg_val /= nodes.size();
        vals[0] = avg_val;
        MoFEMFunctionReturnHot(0);
      };
    } else {
      get_data_to_node = [this, th, &vals](EntityHandle tet) {
        MoFEMFunctionBeginHot;
        CHKERR mwlsMoab.tag_get_data(th, &tet, 1, &*vals.begin());
        MoFEMFunctionReturnHot(0);
      };
    }
 
    // clean tag on nodes
    CHKERR mwlsMoab.tag_clear_data(th, levelNodes, &*vals.begin());
    for (auto tet : levelTets) {
 
      CHKERR get_data_to_node(tet);
      int num_nodes;
      const EntityHandle *conn;
      CHKERR mwlsMoab.get_connectivity(tet, conn, num_nodes, false);
      EntityHandle ho_node;
      EntityType tet_type = mwlsMoab.type_from_handle(tet);
      CHKERR mwlsMoab.high_order_node(tet, conn, tet_type, ho_node);
      CHKERR mwlsMoab.tag_set_data(th, &ho_node, 1, &*vals.begin());
      CHKERR mwlsMoab.tag_set_data(thMidNode, &tet, 1, &ho_node);
    }
 
  } else {
 
    PetscBool add_analytical_internal_stress_operators = PETSC_FALSE;
    CHKERR PetscOptionsGetBool(PETSC_NULL, "",
                               "-add_analytical_internal_stress_operators",
                               &add_analytical_internal_stress_operators, NULL);
 
    CHKERR mwlsMoab.tag_clear_data(th, levelNodes, &*vals.begin());
    for (auto tet : levelTets) {
      int num_nodes;
      const EntityHandle *conn;
      CHKERR mwlsMoab.get_connectivity(tet, conn, num_nodes, false);
      EntityHandle ho_node;
      EntityType tet_type = mwlsMoab.type_from_handle(tet);
      CHKERR mwlsMoab.high_order_node(tet, conn, tet_type, ho_node);
      CHKERR mwlsMoab.tag_set_data(thMidNode, &tet, 1, &ho_node);
 
      if (add_analytical_internal_stress_operators == PETSC_FALSE) {
        VectorDouble coords(3);
        CHKERR mwlsMoab.get_coords(&ho_node, 1, &coords[0]);
        FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3> t_coords(
            &coords[0], &coords[1], &coords[2]);
        auto t_strain = CrackPropagation::analyticalStrainFunction(t_coords);
 
        constexpr double K = 166.666666666e9; // K = E / 3 / (1 - 2*nu),
                                              // E = 200 GPa, nu = 0.3
 
        // Only diagonal, vectorial notation
        vals[0] = 3 * K * t_strain(0, 0);
        vals[1] = 3 * K * t_strain(1, 1);
        vals[2] = 3 * K * t_strain(2, 2);
      }
 
      CHKERR mwlsMoab.tag_set_data(th, &ho_node, 1, &*vals.begin());
    }
  }
 
  MoFEMFunctionReturn(0);
}

◆ loadMWLSMesh()

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

Load mesh with data and do basic calculation

Parameters

file_name File name

Returns: Error code

Definition at line 139 of file MWLS.cpp.

                                                           {
  MoFEMFunctionBegin;
 
  CHKERR getMWLSOptions();
 
  const char *option;
  option = "";
  CHKERR mwlsMoab.load_file(file_name.c_str(), 0, option);
  CHKERR mwlsMoab.get_entities_by_dimension(0, 3, levelTets);
  if (levelTets.empty()) {
    SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
            "No 3d element in MWLS mesh");
  }
 
  CHKERR mwlsMoab.get_adjacencies(levelTets, 1, true, levelEdges,
                                  moab::Interface::UNION);
 
  // Create HO node on tetrahedral
  if (!useNodalData) {
    EntityHandle meshset;
    CHKERR mwlsMoab.create_meshset(MESHSET_SET, meshset);
    CHKERR mwlsMoab.add_entities(meshset, levelTets);
    CHKERR mwlsMoab.convert_entities(meshset, false, false, true);
    CHKERR mwlsMoab.delete_entities(&meshset, 1);
 
    // Get HO nodes in TET center
    std::vector<double> edge_length;
    for (auto tet : levelTets) {
      int num_nodes;
      const EntityHandle *conn;
      CHKERR mwlsMoab.get_connectivity(tet, conn, num_nodes, false);
      EntityHandle ho_node;
      EntityType tit_type = mwlsMoab.type_from_handle(tet);
      CHKERR mwlsMoab.high_order_node(tet, conn, tit_type, ho_node);
      levelNodes.insert(ho_node);
    }
 
  } else {
    Range tets;
    CHKERR mwlsMoab.get_entities_by_dimension(0, 3, tets);
    CHKERR mwlsMoab.get_connectivity(tets, levelNodes, true);
  }
 
  // Store edge nodes coordinates in FTensor
  double edge_node_coords[6];
  FTensor::Tensor1<double *, 3> t_node_edge[2] = {
      FTensor::Tensor1<double *, 3>(edge_node_coords, &edge_node_coords[1],
                                    &edge_node_coords[2]),
      FTensor::Tensor1<double *, 3>(&edge_node_coords[3], &edge_node_coords[4],
                                    &edge_node_coords[5])};
 
  FTensor::Index<'i', 3> i;
 
  // Get edge lengths
  maxEdgeL = 0;
  for (auto edge : levelEdges) {
    int num_nodes;
    const EntityHandle *conn;
    CHKERR mwlsMoab.get_connectivity(edge, conn, num_nodes, true);
    CHKERR mwlsMoab.get_coords(conn, num_nodes, edge_node_coords);
    t_node_edge[0](i) -= t_node_edge[1](i);
    double l = sqrt(t_node_edge[0](i) * t_node_edge[0](i));
    maxEdgeL = (maxEdgeL < l) ? l : maxEdgeL;
  }
 
  // Create tree and find maximal edge length. Maximal edge length is used
  // to find neighbours nodes to material point.
  myTree = boost::make_shared<BVHTree>(&mwlsMoab);
  {
    std::ostringstream opts;
    opts << "MAX_DEPTH=" << maxThreeDepth
        //  << ";MAX_PER_LEAF=" << maxElemsInDMFactor * nbBasePolynomials + 1
         << ";SPLITS_PER_DIR=" << splitsPerDir;
    FileOptions tree_opts(opts.str().c_str());
 
    MOFEM_LOG_CHANNEL("MWLSWorld");
    MOFEM_LOG_TAG("MWLSWorld", "MWLS");
    BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
    MOFEM_LOG("MWLSWorld", Sev::verbose) << "Build tree ... ";
    CHKERR myTree->build_tree(levelTets, &treeRoot, &tree_opts);
    MOFEM_LOG("MWLSWorld", Sev::verbose) << "done";
  }
 
  // moab::BoundBox box;
  // CHKERR myTree->get_bounding_box(box);
  // cout << box << endl;
  // myTree->print();
 
  auto map_elems_positions = [&]() {
    std::map<EntityHandle, std::array<double, 3>> elem_coords_map;
    Range tets;
    CHKERR mwlsMoab.get_entities_by_dimension(0, 3, tets, false);
    MatrixDouble elem_coords;
    for (auto p = tets.pair_begin(); p != tets.pair_end(); ++p) {
      EntityHandle *conn;
      int count, verts_per_e;
      Range slot(p->first, p->second);
      CHKERR mwlsMoab.connect_iterate(slot.begin(), slot.end(), conn,
                                      verts_per_e, count);
      if (count != (int)slot.size())
        THROW_MESSAGE("Should be one continuos sequence");
      elem_coords.resize(verts_per_e, 3, false);
 
      for (auto t = 0; t != count; ++t) {
        CHKERR mwlsMoab.get_coords(conn, verts_per_e,
                                   &*elem_coords.data().begin());
        auto get_center = [&]() {
          std::array<double, 3> c{0, 0, 0};
          for (auto d : {0, 1, 2})
            for (auto n = 0; n < verts_per_e; ++n)
              c[d] += elem_coords(n, d);
          for (auto d : {0, 1, 2})
            c[d] /= verts_per_e;
          return c;
        };
        elem_coords_map.emplace(std::make_pair(p->first + t, get_center()));
        conn += verts_per_e;
      }
    }
    return elem_coords_map;
  };
 
  elemCentreMap = map_elems_positions();
 
  // Set tag for MWLS stress when SSLV116 test is run
  if (mwlsAnalyticalInternalStressTest) {
    Tag th;
    std::array<double, 9> def;
    def.fill(0);
    CHKERR mwlsMoab.tag_get_handle("MED_SSLV116", 9, MB_TYPE_DOUBLE, th,
                                   MB_TAG_CREAT | MB_TAG_SPARSE, def.data());
  }
 
  MoFEMFunctionReturn(0);
}

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:

users_modules/fracture_mechanics/src/MWLS.hpp
users_modules/fracture_mechanics/src/impl/MWLS.cpp

Classes

Public Types

Public Member Functions

Public Attributes

Private Member Functions

Private Attributes

Detailed Description

Member Typedef Documentation

◆ MatDiffDiffVals

◆ MatDiffVals

◆ OpMWLSCalculateBaseCoeffcientsAtGaussPts

◆ OpMWLSCalculateBaseCoeffcientsForFaceAtGaussPts

◆ VecVals

Constructor & Destructor Documentation

◆ MWLSApprox()

◆ ~MWLSApprox()

Member Function Documentation

◆ calculateApproxCoeff()

◆ calculateBase()

◆ calculateDiffBase()

◆ calculateDiffDiffBase()

◆ evalDiffDiffWeight()

◆ evalDiffWeight()

◆ evaluateApproxFun()

◆ evaluateFirstDiffApproxFun()

◆ evaluateSecondDiffApproxFun()

◆ evalWeight()

◆ getDataApprox()

◆ getDiffDataApprox()

◆ getDiffDiffDataApprox()

◆ getErrorAtNode()

◆ getInfluenceNodes() [1/2]

◆ getInfluenceNodes() [2/2]

◆ getInvAB()

◆ getMWLSOptions()

◆ getShortenDm()

◆ getTagData()

◆ getUseGlobalBaseAtMaterialReferenceConfiguration()

◆ getValuesToNodes()

◆ loadMWLSMesh()

Member Data Documentation

◆ A

◆ approxBasePoint

◆ approxFun

◆ approxPointCoords

◆ arcLengthDof

◆ B

◆ baseFun

◆ baseFunInvA

◆ diffA

◆ diffApproxFun

◆ diffB

◆ diffBaseFun

◆ diffDiffApproxFun

◆ diffDiffBaseFun

◆ diffDiffRhoAtGaussPts

◆ diffInvA

◆ diffRhoAtGaussPts

◆ diffStressAtGaussPts

◆ distLeafs

◆ dmFactor

◆ dmNodesMap

◆ elemCentreMap

◆ F_lambda

◆ functionTestHack

◆ influenceNodes

◆ influenceNodesData

◆ influenceNodesMap

◆ invA

◆ invAB

◆ invABMap

◆ L

◆ leafsDist

◆ leafsTetsCentre

◆ leafsTetsVecLeaf

◆ leafsVec

◆ levelEdges

◆ levelNodes

◆ levelTets

◆ maxEdgeL