|
| v0.14.0
|
#include <users_modules/fracture_mechanics/src/MWLS.hpp>
|
| 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 () |
|
|
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 |
|
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.
◆ MatDiffDiffVals
◆ MatDiffVals
◆ OpMWLSCalculateBaseCoeffcientsAtGaussPts
◆ OpMWLSCalculateBaseCoeffcientsForFaceAtGaussPts
◆ VecVals
◆ 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.
45 if (!LogManager::checkIfChannelExist(
"MWLSWorld")) {
46 auto core_log = logging::core::get();
49 LogManager::createSink(LogManager::getStrmWorld(),
"MWLSWorld"));
51 LogManager::createSink(LogManager::getStrmSync(),
"MWLSSync"));
53 LogManager::createSink(LogManager::getStrmSelf(),
"MWLSSelf"));
55 LogManager::setLog(
"MWLSWorld");
56 LogManager::setLog(
"MWLSSync");
57 LogManager::setLog(
"MWLSSelf");
64 MOFEM_LOG(
"MWLSWorld", Sev::noisy) <<
"MWLS created";
◆ ~MWLSApprox()
virtual FractureMechanics::MWLSApprox::~MWLSApprox |
( |
| ) |
|
|
virtualdefault |
◆ 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.
580 auto to_range = [&](
auto ents) {
582 r.insert_list(ents.begin(), ents.end());
586 auto save_entities = [&](
const std::string name,
Range ents) {
600 auto trim_nodes = [&](
const double dm) {
608 &*node_coors_vec.begin());
610 &node_coors_vec[0], &node_coors_vec[1], &node_coors_vec[2]);
614 t_node_coords(
i) -= t_approx_point(
i);
615 const double r2 = t_node_coords(
i) * t_node_coords(
i);
617 influence_nodes_tmp.emplace_back(node);
624 auto eval_AB = [
this, global_derivatives](
const double dm) {
632 for (
int d = 0;
d != 3;
d++) {
644 &*nodes_coords.begin());
645 double *ptr = &*nodes_coords.begin();
647 ptr, ptr + 1, ptr + 2);
649 double min_distance = -1;
653 t_node_coords(
i) -= t_approx_point(
i);
655 const double r = sqrt(t_node_coords(
i) * t_node_coords(
i));
659 if (!ii ||
r < min_distance) {
666 t_diff_r(
i) = (1. /
r) * t_node_coords(
i) * (-1 /
shortenDm);
672 calculateBase<FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>>(
686 for (
int j = 0;
j <=
k;
j++) {
690 if (global_derivatives) {
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++) {
698 diffB[
d](
k, ii) = diff_wp_r_dx;
706 "Nearest node not found \n");
713 auto invert_A = [
this]() {
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);
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);
735 for (
size_t i = 0;
i != nb_base_functions; ++
i)
736 for (
size_t j = 0;
j != nb_base_functions; ++
j)
744 auto throw_error = [&]() {
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 "
761 auto solve_one = [&]() {
783 auto calculate_shape_functions_coefficients = [
this]() {
796 CHKERR calculate_shape_functions_coefficients();
◆ 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.
689 baseFun[4] = t_coords(0) * t_coords(0);
690 baseFun[5] = t_coords(1) * t_coords(1);
691 baseFun[6] = t_coords(2) * t_coords(2);
692 baseFun[7] = t_coords(0) * t_coords(1);
693 baseFun[8] = t_coords(0) * t_coords(2);
694 baseFun[9] = t_coords(1) * t_coords(2);
◆ 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.
705 for (
int d = 0;
d != 3; ++
d) {
710 for (
int d = 0;
d != 3;
d++) {
◆ 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.
754 for (
int d = 0;
d != 9; ++
d) {
763 const double dm2 = dm * dm;
◆ evalDiffDiffWeight()
double FractureMechanics::MWLSApprox::evalDiffDiffWeight |
( |
const double |
r | ) |
|
|
inlineprivate |
Definition at line 812 of file MWLS.hpp.
813 const double rr =
r *
r;
815 return -36 * rr + 48 *
r - 12;
◆ evalDiffWeight()
double FractureMechanics::MWLSApprox::evalDiffWeight |
( |
const double |
r | ) |
|
|
inlineprivate |
Definition at line 797 of file MWLS.hpp.
798 const double rr =
r *
r;
799 const double rrr = rr *
r;
801 return -12 *
r + 24 * rr - 12 * rrr;
◆ evaluateApproxFun()
evaluate approximation function at material point
- Parameters
-
material_coords | material position |
- Returns
- error code
Definition at line 801 of file MWLS.cpp.
805 auto cal_base_functions_at_eval_point = [
this, eval_point](
const double dm) {
815 t_delta(
i) = t_eval_point(
i) - t_approx_point(
i);
822 auto calculate_shape_functions = [
this]() {
831 CHKERR calculate_shape_functions();
◆ 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.
840 auto cal_diff_base_functions_at_eval_point = [
this,
841 eval_point](
const double dm) {
851 t_delta(
i) = t_eval_point(
i) - t_approx_point(
i);
858 auto calculate_global_shape_functions_derivatives = [
this]() {
867 for (
int d = 0;
d != 3;
d++) {
872 for (
int d = 0;
d != 3;
d++) {
878 for (
int d = 0;
d != 3; ++
d) {
896 auto calculate_local_shape_functions_derivatives = [
this]() {
902 for (
int d = 0;
d != 3; ++
d) {
912 if (global_derivatives)
913 CHKERR calculate_global_shape_functions_derivatives();
915 CHKERR calculate_local_shape_functions_derivatives();
◆ 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.
925 auto cal_diff_diff_base_functions_at_eval_point =
926 [
this, eval_point](
const double dm) {
936 t_delta(
i) = t_eval_point(
i) - t_approx_point(
i);
943 auto calculate_global_shape_functions_second_derivatives = [
this]() {
949 auto calculate_local_shape_functions_second_derivatives = [
this]() {
957 for (
int n = 0;
n != 3; ++
n) {
958 for (
int d = 0;
d != 3; ++
d) {
959 const int indx =
d + 3 *
n;
969 if (global_derivatives) {
973 CHKERR calculate_local_shape_functions_second_derivatives();
◆ evalWeight()
double FractureMechanics::MWLSApprox::evalWeight |
( |
const double |
r | ) |
|
|
inlineprivate |
Definition at line 781 of file MWLS.hpp.
782 const double rr =
r *
r;
783 const double rrr = rr *
r;
784 const double r4 = rrr *
r;
786 return 1 - 6 * rr + 8 * rrr - 3 * r4;
◆ getDataApprox()
Get values at point
- Returns
- error code
Definition at line 992 of file MWLS.cpp.
◆ getDiffDataApprox()
Get derivatives of values at points
- Returns
- error code
Definition at line 997 of file MWLS.cpp.
◆ getDiffDiffDataApprox()
Get second derivatives of values at points
- Returns
- error code
Definition at line 1002 of file MWLS.cpp.
◆ 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.
1022 auto get_compliance_matrix = [&]() {
1024 t_C(
i,
j,
k,
l) = 0;
1028 t_C(0, 0, 1, 1) = -
c;
1029 t_C(0, 0, 2, 2) = -
c;
1031 t_C(1, 1, 0, 0) = -
c;
1033 t_C(1, 1, 2, 2) = -
c;
1035 t_C(2, 2, 0, 0) = -
c;
1036 t_C(2, 2, 1, 1) = -
c;
1040 t_C(0, 1, 0, 1) =
d;
1041 t_C(0, 2, 0, 2) =
d;
1042 t_C(1, 2, 1, 2) =
d;
1047 auto t_C = get_compliance_matrix();
1051 "Tag length is not consistent with outData size. \nMost probably "
1052 "getTagData has not been invoked yet!\n");
1055 val_at_nearest_influence_node.resize(length,
false);
1056 val_at_nearest_influence_node.clear();
1058 &*val_at_nearest_influence_node.begin());
1062 auto get_symm_tensor = [](
auto &
m) {
1065 &
m(0), &
m(3), &
m(4),
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);
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));
1083 auto get_stress_nrm2 = [
i,
j](
auto &t_s) {
return t_s(
i,
j) * t_s(
i,
j); };
1085 total_stress_at_node = get_stress_nrm2(t_stress_at_nearset_node);
1086 total_energy = get_energy_norm(t_stress_at_nearset_node);
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);
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;
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);
◆ getInfluenceNodes() [1/2]
std::vector<EntityHandle>& FractureMechanics::MWLSApprox::getInfluenceNodes |
( |
| ) |
|
|
inline |
◆ getInfluenceNodes() [2/2]
Get node in influence domain.
- Parameters
-
material_coords | material position |
- Returns
- error code
Definition at line 373 of file MWLS.cpp.
374 constexpr
double eps_overlap = 0.01;
378 using DistEntPair = std::pair<double, EntityHandle>;
379 using DistEntMap = multi_index_container<
383 ordered_non_unique<member<DistEntPair, double, &DistEntPair::first>>
395 for (
int dd : {0, 1, 2})
401 auto find_leafs = [&](
const auto dm) {
413 [](
const auto &
a,
const auto &b) { return a.first < b.first; });
417 auto get_influence_nodes = [&](
const auto &tets) {
418 std::vector<EntityHandle> influence_nodes;
420 influence_nodes.resize(tets.size());
422 &*influence_nodes.begin());
425 CHKERR mwlsMoab.get_connectivity(&*tets.begin(), tets.size(), nodes,
427 influence_nodes.clear();
428 influence_nodes.reserve(nodes.size());
430 influence_nodes.emplace_back(
n);
432 return influence_nodes;
435 auto to_range = [&](
auto ents) {
437 r.insert_list(ents.begin(), ents.end());
441 auto save_entities = [&](
const std::string name,
Range ents) {
457 std::vector<EntityHandle> leafs_tets;
460 auto get_dm2_from_influence_points = [&]() {
465 &*node_coors_vec.begin());
467 material_coords[0], material_coords[1], material_coords[2]);
469 &node_coors_vec[0], &node_coors_vec[1], &node_coors_vec[2]);
472 t_node_coords(
i) -= t_approx_point(
i);
473 const double r2 = t_node_coords(
i) * t_node_coords(
i);
482 double set_dm2 = shorten_dm2;
483 auto min_max_dist_it = dist_map.get<0>().begin();
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()))) {
503 for (
auto ii : {0, 1, 2}) {
504 t_c(ii) =
c[ii] - material_coords[ii];
517 if (dist_map.size() > max_nb_elems)
518 set_dm2 = std::min(min_max_dist_it->first, shorten_dm2);
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);
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(),
528 for (; nb < max_nb_elems; ++nb)
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)
550 for (
auto &it : dist_map)
552 <<
"dm map " << it.first <<
" " << it.second;
554 MOFEM_LOG(
"MWLSWorld", Sev::error) <<
"leafs found " << leafs.size();
556 <<
"dist map size " << dist_map.size();
557 MOFEM_LOG(
"MWLSWorld", Sev::error) <<
"dm2 " << dm2;
559 MOFEM_LOG_C(
"MWLSWorld", Sev::error,
"point: ( %g %g %g )",
566 const auto dm2_from_influence_pts = get_dm2_from_influence_points();
567 shortenDm = (1 + eps_overlap) * sqrt(dm2_from_influence_pts);
◆ getInvAB()
MatrixDouble& FractureMechanics::MWLSApprox::getInvAB |
( |
| ) |
|
|
inline |
◆ getMWLSOptions()
get options from command line for MWLS
- Returns
- error code
Definition at line 71 of file MWLS.cpp.
74 "Get options for moving least square approximation",
77 CHKERR PetscOptionsScalar(
"-mwls_dm",
78 "edge length scaling for influence radius",
"",
80 MOFEM_LOG_C(
"MWLSWorld", Sev::inform,
"### Input parameter: -mwls_dm %6.4e",
84 "-mwls_number_of_base_functions",
85 "1 = constant, 4 = linear, 10 = quadratic approximation",
"",
88 "### Input parameter: -mwls_number_of_base_functions %d",
91 CHKERR PetscOptionsInt(
"-mwls_max_elems_factor",
92 "Set max elements factor max_nb_elems = "
93 "maxElemsInDMFactor * nbBasePolynomials",
97 "### Input parameter: -mwls_max_elems_factor %d",
101 "-mwls_use_global_base_at_reference_configuration",
102 "true local mwls base functions at reference configuration are used",
"",
108 "-mwls_use_nodal_data",
109 "true nodal data values are used (without averaging to the midpoint)",
"",
112 "### Input parameter: -mwls_use_nodal_data %d",
useNodalData);
115 "-mwls_analytical_internal_stress",
116 "use analytical strain field for internal stress mwls test",
"",
120 "### Input parameter: -test_mwls_internal_stress_analytical %d",
124 CHKERR PetscOptionsInt(
"-mwls_tree_depth",
"maximal three depths",
"",
129 CHKERR PetscOptionsInt(
"-mwls_splits_per_dir",
"splits per direction",
"",
134 ierr = PetscOptionsEnd();
◆ getShortenDm()
auto& FractureMechanics::MWLSApprox::getShortenDm |
( |
| ) |
|
|
inline |
◆ getTagData()
Get tag data on influence nodes
- Parameters
-
- Returns
- error code
Definition at line 978 of file MWLS.cpp.
◆ 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.
◆ 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.
281 MB_TAG_CREAT | MB_TAG_DENSE);
295 if (name.compare(
"RHO") == 0) {
301 for (
auto &
n : nodes) {
306 avg_val /= nodes.size();
322 CHKERR get_data_to_node(tet);
327 EntityType tet_type =
mwlsMoab.type_from_handle(tet);
335 PetscBool add_analytical_internal_stress_operators = PETSC_FALSE;
337 "-add_analytical_internal_stress_operators",
338 &add_analytical_internal_stress_operators, NULL);
346 EntityType tet_type =
mwlsMoab.type_from_handle(tet);
350 if (add_analytical_internal_stress_operators == PETSC_FALSE) {
354 &coords[0], &coords[1], &coords[2]);
357 constexpr
double K = 166.666666666e9;
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);
◆ loadMWLSMesh()
MoFEMErrorCode FractureMechanics::MWLSApprox::loadMWLSMesh |
( |
std::string |
file_name | ) |
|
Load mesh with data and do basic calculation
- Parameters
-
- Returns
- Error code
Definition at line 139 of file MWLS.cpp.
150 "No 3d element in MWLS mesh");
154 moab::Interface::UNION);
165 std::vector<double> edge_length;
171 EntityType tit_type =
mwlsMoab.type_from_handle(tet);
183 double edge_node_coords[6];
186 &edge_node_coords[2]),
188 &edge_node_coords[5])};
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));
208 std::ostringstream opts;
212 FileOptions tree_opts(opts.str().c_str());
216 BOOST_LOG_SCOPED_THREAD_ATTR(
"Timeline", attrs::timer());
217 MOFEM_LOG(
"MWLSWorld", Sev::verbose) <<
"Build tree ... ";
219 MOFEM_LOG(
"MWLSWorld", Sev::verbose) <<
"done";
227 auto map_elems_positions = [&]() {
228 std::map<EntityHandle, std::array<double, 3>> elem_coords_map;
232 for (
auto p = tets.pair_begin(); p != tets.pair_end(); ++p) {
234 int count, verts_per_e;
235 Range slot(p->first, p->second);
238 if (count != (
int)slot.size())
240 elem_coords.resize(verts_per_e, 3,
false);
242 for (
auto t = 0;
t != count; ++
t) {
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})
254 elem_coords_map.emplace(std::make_pair(p->first +
t, get_center()));
258 return elem_coords_map;
266 std::array<double, 9> def;
269 MB_TAG_CREAT | MB_TAG_SPARSE, def.data());
ublas::symmetric_matrix<double, ublas::lower> FractureMechanics::MWLSApprox::A |
|
private |
◆ approxBasePoint
MatrixDouble FractureMechanics::MWLSApprox::approxBasePoint |
◆ approxFun
VectorDouble FractureMechanics::MWLSApprox::approxFun |
|
private |
◆ approxPointCoords
VectorDouble3 FractureMechanics::MWLSApprox::approxPointCoords |
◆ arcLengthDof
boost::weak_ptr<DofEntity> FractureMechanics::MWLSApprox::arcLengthDof |
MatrixDouble FractureMechanics::MWLSApprox::B |
|
private |
◆ baseFun
VectorDouble FractureMechanics::MWLSApprox::baseFun |
|
private |
◆ baseFunInvA
VectorDouble FractureMechanics::MWLSApprox::baseFunInvA |
|
private |
◆ diffA
ublas::symmetric_matrix<double, ublas::lower, ublas::row_major, DoubleAllocator> FractureMechanics::MWLSApprox::diffA[3] |
|
private |
◆ diffApproxFun
MatrixDouble FractureMechanics::MWLSApprox::diffApproxFun |
|
private |
◆ diffB
MatrixDouble FractureMechanics::MWLSApprox::diffB[3] |
|
private |
◆ diffBaseFun
VectorDouble FractureMechanics::MWLSApprox::diffBaseFun[3] |
|
private |
◆ diffDiffApproxFun
MatrixDouble FractureMechanics::MWLSApprox::diffDiffApproxFun |
|
private |
◆ diffDiffBaseFun
VectorDouble FractureMechanics::MWLSApprox::diffDiffBaseFun[9] |
|
private |
◆ diffDiffRhoAtGaussPts
boost::shared_ptr<MatrixDouble> FractureMechanics::MWLSApprox::diffDiffRhoAtGaussPts |
◆ diffInvA
ublas::symmetric_matrix<double, ublas::lower, ublas::row_major, DoubleAllocator> FractureMechanics::MWLSApprox::diffInvA[3] |
|
private |
◆ diffRhoAtGaussPts
boost::shared_ptr<MatrixDouble> FractureMechanics::MWLSApprox::diffRhoAtGaussPts |
◆ diffStressAtGaussPts
MatrixDouble FractureMechanics::MWLSApprox::diffStressAtGaussPts |
◆ distLeafs
◆ dmFactor
double FractureMechanics::MWLSApprox::dmFactor |
|
private |
Relative value of weight function radius.
Definition at line 617 of file MWLS.hpp.
◆ 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 |
◆ F_lambda
Vec FractureMechanics::MWLSApprox::F_lambda |
◆ 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 |
◆ influenceNodesData
MatrixDouble FractureMechanics::MWLSApprox::influenceNodesData |
|
private |
◆ influenceNodesMap
Influence nodes on elements on at integration points.
Definition at line 230 of file MWLS.hpp.
◆ invA
MatrixDouble FractureMechanics::MWLSApprox::invA |
|
private |
◆ invAB
MatrixDouble FractureMechanics::MWLSApprox::invAB |
|
private |
◆ 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.
ublas::triangular_matrix<double, ublas::lower> FractureMechanics::MWLSApprox::L |
|
private |
◆ leafsDist
std::vector<double> FractureMechanics::MWLSApprox::leafsDist |
|
private |
◆ leafsTetsCentre
std::vector<double> FractureMechanics::MWLSApprox::leafsTetsCentre |
|
private |
◆ leafsTetsVecLeaf
std::vector<EntityHandle> FractureMechanics::MWLSApprox::leafsTetsVecLeaf |
|
private |
◆ leafsVec
std::vector<EntityHandle> FractureMechanics::MWLSApprox::leafsVec |
|
private |
◆ levelEdges
Range FractureMechanics::MWLSApprox::levelEdges |
|
private |
◆ levelNodes
Range FractureMechanics::MWLSApprox::levelNodes |
|
private |
◆ levelTets
Range FractureMechanics::MWLSApprox::levelTets |
|
private |
◆ maxEdgeL
double FractureMechanics::MWLSApprox::maxEdgeL |
|
private |
◆ 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
◆ mwlsAnalyticalInternalStressTest
PetscBool FractureMechanics::MWLSApprox::mwlsAnalyticalInternalStressTest |
|
private |
◆ mwlsCore
moab::Core FractureMechanics::MWLSApprox::mwlsCore |
◆ mwlsMaterialPositions
MatrixDouble FractureMechanics::MWLSApprox::mwlsMaterialPositions |
◆ mwlsMoab
moab::Interface& FractureMechanics::MWLSApprox::mwlsMoab |
◆ myTree
boost::shared_ptr<BVHTree> FractureMechanics::MWLSApprox::myTree |
|
private |
◆ 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 |
◆ outData
VecVals FractureMechanics::MWLSApprox::outData |
|
private |
◆ outDeviatoricDifference
VecVals FractureMechanics::MWLSApprox::outDeviatoricDifference |
|
private |
◆ outDeviatoricError
double FractureMechanics::MWLSApprox::outDeviatoricError |
|
private |
◆ outDiffData
◆ outDiffDiffData
◆ outHydroStaticError
double FractureMechanics::MWLSApprox::outHydroStaticError |
|
private |
◆ rhoAtGaussPts
boost::shared_ptr<VectorDouble> FractureMechanics::MWLSApprox::rhoAtGaussPts |
◆ 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 |
◆ singularInitialDisplacement
MatrixDouble FractureMechanics::MWLSApprox::singularInitialDisplacement |
◆ 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 |
◆ testA
ublas::symmetric_matrix<double, ublas::lower> FractureMechanics::MWLSApprox::testA |
|
private |
◆ tetsInBlock
Range FractureMechanics::MWLSApprox::tetsInBlock |
◆ thMidNode
Tag FractureMechanics::MWLSApprox::thMidNode |
◆ treeRoot
◆ treeTets
std::vector<EntityHandle> FractureMechanics::MWLSApprox::treeTets |
|
private |
◆ trisInBlock
Range FractureMechanics::MWLSApprox::trisInBlock |
◆ useGlobalBaseAtMaterialReferenceConfiguration
PetscBool FractureMechanics::MWLSApprox::useGlobalBaseAtMaterialReferenceConfiguration |
|
private |
◆ useNodalData
PetscBool FractureMechanics::MWLSApprox::useNodalData |
|
private |
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
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
double evalWeight(const double r)
double evalDiffWeight(const double r)
std::map< EntityHandle, std::array< double, 3 > > elemCentreMap
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
double young_modulus
Young modulus.
MatrixDouble diffDiffApproxFun
MoFEMErrorCode getMWLSOptions()
get options from command line for MWLS
virtual MPI_Comm & get_comm() const =0
boost::shared_ptr< BVHTree > myTree
int splitsPerDir
Split per direction in the tree.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
boost::shared_ptr< MatrixDouble > diffRhoAtGaussPts
UBlasMatrix< double > MatrixDouble
MoFEMErrorCode calculateBase(const T &t_coords)
MatrixDouble diffApproxFun
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
std::vector< EntityHandle > leafsVec
std::vector< EntityHandle > leafsTetsVecLeaf
VectorDouble3 approxPointCoords
ublas::triangular_matrix< double, ublas::lower > L
MatrixDouble influenceNodesData
std::vector< EntityHandle > influenceNodes
const double c
speed of light (cm/ns)
#define CHKERR
Inline error check.
ublas::symmetric_matrix< double, ublas::lower, ublas::row_major, DoubleAllocator > diffInvA[3]
boost::shared_ptr< VectorDouble > rhoAtGaussPts
boost::weak_ptr< DofEntity > arcLengthDof
#define MOFEM_LOG_C(channel, severity, format,...)
ublas::symmetric_matrix< double, ublas::lower > A
std::vector< double > leafsDist
double dmFactor
Relative value of weight function radius.
double poisson_ratio
Poisson ratio.
int maxThreeDepth
Maximal three depths.
moab::Interface & mwlsMoab
ublas::symmetric_matrix< double, ublas::lower > testA
static PetscErrorCode solve(Mat mat, Vec x, Vec y)
@ MOFEM_OPERATION_UNSUCCESSFUL
MoFEMErrorCode calculateDiffDiffBase(const T &t_coords, const double dm)
MatDiffDiffVals outDiffDiffData
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
constexpr double t
plate stiffness
FTensor::Index< 'i', SPACE_DIM > i
static FTensor::Tensor2_symmetric< double, 3 > analyticalStrainFunction(FTensor::Tensor1< FTensor::PackPtr< double *, 1 >, 3 > &t_coords)
PetscBool mwlsAnalyticalInternalStressTest
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)
#define MOFEM_LOG(channel, severity)
Log.
MoFEM::Interface & mField
ublas::symmetric_matrix< double, ublas::lower, ublas::row_major, DoubleAllocator > diffA[3]
std::vector< double > leafsTetsCentre
VectorDouble diffBaseFun[3]
FTensor::Index< 'j', 3 > j
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
MoFEMErrorCode calculateDiffBase(const T &t_coords, const double dm)
@ MOFEM_DATA_INCONSISTENCY
size_t cholesky_decompose(const MATRIX &A, TRIA &L)
decompose the symmetric positive definit matrix A into product L L^T.
void cholesky_solve(const TRIA &L, VEC &x, ublas::lower)
solve system L L^T x = b inplace
UBlasVector< double > VectorDouble
FTensor::Index< 'm', 3 > m
Kronecker Delta class symmetric.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
PetscBool useGlobalBaseAtMaterialReferenceConfiguration
std::vector< EntityHandle > treeTets
FTensor::Index< 'k', 3 > k
EntityHandle nearestInfluenceNode
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
boost::shared_ptr< MatrixDouble > diffDiffRhoAtGaussPts
VectorDouble diffDiffBaseFun[9]
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
FTensor::Index< 'l', 3 > l
std::vector< std::pair< double, EntityHandle > > distLeafs
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)