Skip to content
Snippets Groups Projects
Commit a8adcdeb authored by Dmitri Naumov's avatar Dmitri Naumov
Browse files

Merge branch 'material_anisotropy' into 'master'

To implement a decomposition of strain energy density in phase field approach to fracture accounting for orthotropic materials

See merge request ogs/ogs!4215
parents c2fbe319 233eb2fa
No related branches found
No related tags found
No related merge requests found
Showing
with 859 additions and 60 deletions
......@@ -2,3 +2,5 @@ Implemented energy split models are:
- isotropic model \cite bourdin2000numerical, (energy_split_model="Isotropic")
- volumetric-deviatoric split model \cite amor2009regularized (energy_split_model="VolumetricDeviatoric")
- effective stress model \cite wu2017unified (energy_split_model="EffectiveStress")
- orthotropic volumetric-deviatoric model \cite ziaei2022orthogonal, (energy_split_model="OrthoVolDev")
- orthotropic no-tension model \cite ziaei2022orthogonal, (energy_split_model="OrthoMasonry")
......@@ -431,3 +431,14 @@ URL = {https://doi.org/10.1680/geot.2008.58.3.157}
volume = "19",
year = "1979"
}
@misc{ziaei2022orthogonal,
doi = {10.48550/ARXIV.2207.01557},
url = {https://arxiv.org/abs/2207.01557},
author = {Ziaei-Rad, Vahid and Mollaali, Mostafa and Nagel, Thomas and Kolditz, Olaf and Yoshioka, Keita},
keywords = {Numerical Analysis (math.NA), FOS: Mathematics, FOS: Mathematics},
title = {Orthogonal decomposition of anisotropic constitutive models for the phase field approach to fracture},
publisher = {arXiv},
year = {2022},
copyright = {Creative Commons Attribution Non Commercial No Derivatives 4.0 International}
}
\ No newline at end of file
......@@ -50,10 +50,7 @@ std::tuple<MathLib::KelvinVector::KelvinVectorType<
DisplacementDim> /* sigma_real */,
MathLib::KelvinVector::KelvinVectorType<
DisplacementDim> /* sigma_tensile */,
MathLib::KelvinVector::KelvinMatrixType<
DisplacementDim> /* C_tensile */,
MathLib::KelvinVector::KelvinMatrixType<
DisplacementDim> /* C_compressive */,
MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> /* D */,
double /* strain_energy_tensile */, double /* elastic_energy */
>
calculateVolDevDegradedStress(
......@@ -113,8 +110,11 @@ calculateVolDevDegradedStress(
mu * epsd_curr.transpose() * epsd_curr;
KelvinVector const sigma_real =
degradation * sigma_tensile + sigma_compressive;
return std::make_tuple(sigma_real, sigma_tensile, C_tensile, C_compressive,
strain_energy_tensile, elastic_energy);
KelvinMatrix const D = degradation * C_tensile + C_compressive;
return std::make_tuple(sigma_real, sigma_tensile, D, strain_energy_tensile,
elastic_energy);
}
template <int DisplacementDim>
......@@ -122,10 +122,7 @@ std::tuple<MathLib::KelvinVector::KelvinVectorType<
DisplacementDim> /* sigma_real */,
MathLib::KelvinVector::KelvinVectorType<
DisplacementDim> /* sigma_tensile */,
MathLib::KelvinVector::KelvinMatrixType<
DisplacementDim> /* C_tensile */,
MathLib::KelvinVector::KelvinMatrixType<
DisplacementDim> /* C_compressive */,
MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> /* D */,
double /* strain_energy_tensile */, double /* elastic_energy */
>
calculateIsotropicDegradedStress(
......@@ -162,8 +159,10 @@ calculateIsotropicDegradedStress(
C_tensile.noalias() += 2 * mu * P_dev * KelvinMatrix::Identity();
KelvinVector const sigma_real = degradation * sigma_tensile;
return std::make_tuple(sigma_real, sigma_tensile, C_tensile, C_compressive,
strain_energy_tensile, elastic_energy);
KelvinMatrix const D = degradation * C_tensile + C_compressive;
return std::make_tuple(sigma_real, sigma_tensile, D, strain_energy_tensile,
elastic_energy);
}
template <int DisplacementDim>
......@@ -171,10 +170,7 @@ std::tuple<MathLib::KelvinVector::KelvinVectorType<
DisplacementDim> /* sigma_real */,
MathLib::KelvinVector::KelvinVectorType<
DisplacementDim> /* sigma_tensile */,
MathLib::KelvinVector::KelvinMatrixType<
DisplacementDim> /* C_tensile */,
MathLib::KelvinVector::KelvinMatrixType<
DisplacementDim> /* C_compressive */,
MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> /* D */,
double /* strain_energy_tensile */, double /* elastic_energy */
>
calculateIsotropicDegradedStressWithRankineEnergy(
......@@ -227,8 +223,9 @@ calculateIsotropicDegradedStressWithRankineEnergy(
C_tensile.noalias() += 2 * mu * P_dev * KelvinMatrix::Identity();
KelvinVector const sigma_real = degradation * sigma_tensile;
return std::make_tuple(sigma_real, sigma_tensile, C_tensile, C_compressive,
strain_energy_tensile, elastic_energy);
KelvinMatrix const D = degradation * C_tensile + C_compressive;
return std::make_tuple(sigma_real, sigma_tensile, D, strain_energy_tensile,
elastic_energy);
}
} // namespace Phasefield
......
/**
* \file
* \copyright
* Copyright (c) 2012-2022, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*/
#pragma once
#include "MathLib/KelvinVector.h"
namespace MaterialLib
{
namespace Solids
{
namespace Phasefield
{
template <int DisplacementDim>
std::tuple<MathLib::KelvinVector::KelvinVectorType<
DisplacementDim> /* eps_tensile */,
MathLib::KelvinVector::KelvinVectorType<
DisplacementDim> /* sigma_real */,
MathLib::KelvinVector::KelvinVectorType<
DisplacementDim> /* sigma_tensile */,
MathLib::KelvinVector::KelvinVectorType<
DisplacementDim> /* sigma_compressive */,
MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> /* D */,
double /* strain_energy_tensile */, double /* elastic_energy */
>
calculateOrthoVolDevDegradedStress(
double const degradation,
MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const& eps,
MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> const& C_ortho)
{
static constexpr int KelvinVectorSize =
MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
using KelvinVector =
MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
using KelvinMatrix =
MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>;
using Invariants = MathLib::KelvinVector::Invariants<KelvinVectorSize>;
// calculation of deviatoric parts
auto const& P_dev = Invariants::deviatoric_projection;
KelvinMatrix const C = C_ortho;
// calculation of square root of elasticity tensor
Eigen::SelfAdjointEigenSolver<KelvinMatrix> es(C);
KelvinMatrix const sqrt_C = es.operatorSqrt();
Eigen::SelfAdjointEigenSolver<KelvinMatrix> es_inverse(C.inverse());
KelvinMatrix const sqrt_C_inverse = es_inverse.operatorSqrt();
// strain in a transformed space
KelvinVector const epst = sqrt_C * eps;
double const epst_curr_trace = Invariants::trace(epst);
// projection tensors in transformed space
KelvinMatrix teps_p = KelvinMatrix::Zero();
KelvinMatrix teps_n = KelvinMatrix::Zero();
if (epst_curr_trace >= 0) /* QQQ */
{
teps_p.template topLeftCorner<3, 3>().setConstant(1. / 3);
}
else
{
teps_n.template topLeftCorner<3, 3>().setConstant(1. / 3);
}
teps_p.noalias() += P_dev * KelvinMatrix::Identity();
// strain tensile / compressive
KelvinVector const eps_tensile = sqrt_C_inverse * (teps_p * sqrt_C) * eps;
KelvinVector const eps_compressive =
sqrt_C_inverse * (teps_n * sqrt_C) * eps;
// projection tensors in original space
KelvinMatrix const der_eps_p = sqrt_C_inverse * (teps_p * sqrt_C);
KelvinMatrix const der_eps_n = sqrt_C_inverse * (teps_n * sqrt_C);
// C tensile / compressive
KelvinMatrix const C_tensile = der_eps_p.transpose() * C * der_eps_p;
KelvinMatrix const C_compressive = der_eps_n.transpose() * C * der_eps_n;
// stress tensile / compressive
KelvinVector const sigma_tensile = C_tensile * eps;
KelvinVector const sigma_compressive = C_compressive * eps;
// decomposition of strain energy density
double const strain_energy_tensile =
0.5 * sigma_tensile.adjoint() * eps_tensile;
double const strain_energy_compressive =
0.5 * sigma_compressive.adjoint() * eps_compressive;
double const elastic_energy =
degradation * strain_energy_tensile + strain_energy_compressive;
KelvinVector const sigma_real =
degradation * sigma_tensile + sigma_compressive;
KelvinMatrix const D = degradation * C_tensile + C_compressive;
return std::make_tuple(eps_tensile, sigma_real, sigma_tensile,
sigma_compressive, D, strain_energy_tensile,
elastic_energy);
}
template <int DisplacementDim>
std::tuple<MathLib::KelvinVector::KelvinVectorType<
DisplacementDim> /* eps_tensile */,
MathLib::KelvinVector::KelvinVectorType<
DisplacementDim> /* sigma_real */,
MathLib::KelvinVector::KelvinVectorType<
DisplacementDim> /* sigma_tensile */,
MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> /* D */,
double /* strain_energy_tensile */, double /* elastic_energy */
>
calculateOrthoMasonryDegradedStress(
double const degradation,
MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const& eps,
MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> const& C_ortho)
{
using KelvinVector =
MathLib::KelvinVector::KelvinVectorType<DisplacementDim>;
using KelvinMatrix =
MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>;
KelvinMatrix const C = C_ortho;
// calculation of square root of elasticity tensor
Eigen::SelfAdjointEigenSolver<KelvinMatrix> es(C);
KelvinMatrix const sqrt_C = es.operatorSqrt();
Eigen::SelfAdjointEigenSolver<KelvinMatrix> es_inverse(C.inverse());
KelvinMatrix const sqrt_C_inverse = es_inverse.operatorSqrt();
// strain in a transformed space
KelvinVector const epst = sqrt_C * eps;
// strain tensor in 3D
Eigen::Matrix3d const eps_3D =
MathLib::KelvinVector::kelvinVectorToTensor(epst);
Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> const es_eps_3D(eps_3D);
Eigen::Vector3d const eigen_values_eps_3D = es_eps_3D.eigenvalues().real();
Eigen::Matrix3d const eigen_vectors_eps_3D = es_eps_3D.eigenvectors();
Eigen::Matrix3d const E1_eigenp =
eigen_vectors_eps_3D.col(0) * eigen_vectors_eps_3D.col(0).transpose();
Eigen::Matrix3d const E3_eigenp =
eigen_vectors_eps_3D.col(2) * eigen_vectors_eps_3D.col(2).transpose();
KelvinVector const E1_vec =
MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(E1_eigenp);
KelvinVector const E3_vec =
MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(E3_eigenp);
KelvinMatrix const E1oE1 = E1_vec * E1_vec.transpose();
KelvinMatrix const E3oE3 = E3_vec * E3_vec.transpose();
KelvinMatrix I_p = KelvinMatrix::Zero();
KelvinMatrix I_n = KelvinMatrix::Zero();
KelvinMatrix const I_S = KelvinMatrix::Identity();
if (DisplacementDim == 2)
{
if (std::abs(eigen_values_eps_3D(2) - eigen_values_eps_3D(0)) >
std::numeric_limits<double>::epsilon())
{
I_p = (macaulayTensile(eigen_values_eps_3D(0)) -
macaulayTensile(eigen_values_eps_3D(2))) /
(eigen_values_eps_3D(0) - eigen_values_eps_3D(2)) *
(I_S - (E1oE1 + E3oE3)) +
(heaviside(eigen_values_eps_3D(0)) * E1oE1 +
heaviside(eigen_values_eps_3D(2)) * E3oE3);
I_n = (macaulayCompressive(eigen_values_eps_3D(0)) -
macaulayCompressive(eigen_values_eps_3D(2))) /
(eigen_values_eps_3D(0) - eigen_values_eps_3D(2)) *
(I_S - (E1oE1 + E3oE3)) +
(heaviside(-eigen_values_eps_3D(0)) * E1oE1 +
heaviside(-eigen_values_eps_3D(2)) * E3oE3);
}
else
{
I_p = heaviside(eigen_values_eps_3D(0)) * I_S;
I_n = heaviside(-eigen_values_eps_3D(0)) * I_S;
}
}
// strain tensile / compressive
KelvinVector const eps_tensile = sqrt_C_inverse * (I_p * sqrt_C) * eps;
KelvinVector const eps_compressive = sqrt_C_inverse * (I_n * sqrt_C) * eps;
// C tensile / compressive
KelvinMatrix const C_tensile = sqrt_C * I_p * sqrt_C;
KelvinMatrix const C_compressive = sqrt_C * I_n * sqrt_C;
// stress tensile / compressive
KelvinVector const sigma_tensile = C * eps_tensile;
KelvinVector const sigma_compressive = C * eps_compressive;
// decomposition of strain energy density
double const strain_energy_tensile =
0.5 * sigma_tensile.adjoint() * eps_tensile;
double const strain_energy_compressive =
0.5 * sigma_compressive.adjoint() * eps_compressive;
double const elastic_energy =
degradation * strain_energy_tensile + strain_energy_compressive;
KelvinVector const sigma_real =
degradation * sigma_tensile + sigma_compressive;
KelvinMatrix const D = degradation * C_tensile + C_compressive;
return std::make_tuple(eps_tensile, sigma_real, sigma_tensile, D,
strain_energy_tensile, elastic_energy);
}
} // namespace Phasefield
} // namespace Solids
} // namespace MaterialLib
......@@ -247,6 +247,14 @@ std::unique_ptr<Process> createPhaseFieldProcess(
{
return EnergySplitModel::EffectiveStress;
}
if (energy_split_model_string == "OrthoVolDev")
{
return EnergySplitModel::OrthoVolDev;
}
if (energy_split_model_string == "OrthoMasonry")
{
return EnergySplitModel::OrthoMasonry;
}
OGS_FATAL(
"energy_split_model must be 'Isotropic' or 'VolumetricDeviatoric' "
"but '{:s}' was given",
......
......@@ -29,12 +29,30 @@ struct PhaseFieldLocalAssemblerInterface
std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtSigmaTensile(
const double t,
std::vector<GlobalVector*> const& x,
std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtSigmaCompressive(
const double t,
std::vector<GlobalVector*> const& x,
std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtEpsilon(
const double t,
std::vector<GlobalVector*> const& x,
std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
std::vector<double>& cache) const = 0;
virtual std::vector<double> const& getIntPtEpsilonTensile(
const double t,
std::vector<GlobalVector*> const& x,
std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
std::vector<double>& cache) const = 0;
virtual void computeCrackIntegral(
std::size_t mesh_item_id,
std::vector<
......
......@@ -97,8 +97,7 @@ void PhaseFieldLocalAssembler<ShapeFunction, DisplacementDim>::
_process_data.energy_split_model);
auto& sigma = _ip_data[ip].sigma;
auto const& C_tensile = _ip_data[ip].C_tensile;
auto const& C_compressive = _ip_data[ip].C_compressive;
auto const& D = _ip_data[ip].D;
typename ShapeMatricesType::template MatrixType<DisplacementDim,
displacement_size>
......@@ -120,8 +119,7 @@ void PhaseFieldLocalAssembler<ShapeFunction, DisplacementDim>::
(B.transpose() * sigma - N_u.transpose() * rho_sr * b -
local_pressure * N_u.transpose() * dNdx * d) *
w;
local_Jac.noalias() +=
B.transpose() * (degradation * C_tensile + C_compressive) * B * w;
local_Jac.noalias() += B.transpose() * D * B * w;
}
}
......
......@@ -16,12 +16,15 @@
#include "LocalAssemblerInterface.h"
#include "MaterialLib/SolidModels/LinearElasticIsotropic.h"
#include "MaterialLib/SolidModels/LinearElasticIsotropicPhaseField.h"
#include "MaterialLib/SolidModels/LinearElasticOrthotropic.h"
#include "MaterialLib/SolidModels/LinearElasticOrthotropicPhaseField.h"
#include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
#include "NumLib/Fem/InitShapeMatrices.h"
#include "NumLib/Fem/Integration/GenericIntegrationMethod.h"
#include "NumLib/Fem/ShapeMatrixPolicy.h"
#include "ParameterLib/CoordinateSystem.h"
#include "ParameterLib/SpatialPosition.h"
#include "PhaseFieldProcessData.h"
#include "ProcessLib/Deformation/BMatrixPolicy.h"
......@@ -47,7 +50,7 @@ struct IntegrationPointData final
typename ShapeMatrixType::NodalRowVectorType N;
typename ShapeMatrixType::GlobalDimNodalMatrixType dNdx;
typename BMatricesType::KelvinVectorType eps, eps_prev;
typename BMatricesType::KelvinVectorType eps, eps_prev, eps_tensile;
typename BMatricesType::KelvinVectorType sigma_tensile, sigma_compressive,
sigma;
......@@ -58,7 +61,8 @@ struct IntegrationPointData final
DisplacementDim>::MaterialStateVariables>
material_state_variables;
typename BMatricesType::KelvinMatrixType C_tensile, C_compressive;
typename BMatricesType::KelvinMatrixType D;
double integration_weight;
double history_variable, history_variable_prev;
......@@ -90,31 +94,60 @@ struct IntegrationPointData final
{
case EnergySplitModel::Isotropic:
{
std::tie(sigma, sigma_tensile, C_tensile, C_compressive,
strain_energy_tensile, elastic_energy) =
MaterialLib::Solids::Phasefield::
calculateIsotropicDegradedStress<DisplacementDim>(
degradation, bulk_modulus, mu, eps);
std::tie(sigma, sigma_tensile, D, strain_energy_tensile,
elastic_energy) = MaterialLib::Solids::Phasefield::
calculateIsotropicDegradedStress<DisplacementDim>(
degradation, bulk_modulus, mu, eps);
break;
}
case EnergySplitModel::VolDev:
{
std::tie(sigma, sigma_tensile, C_tensile, C_compressive,
strain_energy_tensile, elastic_energy) =
MaterialLib::Solids::Phasefield::
calculateVolDevDegradedStress<DisplacementDim>(
degradation, bulk_modulus, mu, eps);
std::tie(sigma, sigma_tensile, D, strain_energy_tensile,
elastic_energy) = MaterialLib::Solids::Phasefield::
calculateVolDevDegradedStress<DisplacementDim>(
degradation, bulk_modulus, mu, eps);
break;
}
case EnergySplitModel::EffectiveStress:
{
std::tie(sigma, sigma_tensile, C_tensile, C_compressive,
strain_energy_tensile,
std::tie(sigma, sigma_tensile, D, strain_energy_tensile,
elastic_energy) = MaterialLib::Solids::Phasefield::
calculateIsotropicDegradedStressWithRankineEnergy<
DisplacementDim>(degradation, bulk_modulus, mu, eps);
break;
}
case EnergySplitModel::OrthoVolDev:
{
double temperature = 0.;
MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>
C_ortho = static_cast<
MaterialLib::Solids::LinearElasticOrthotropic<
DisplacementDim> const&>(solid_material)
.getElasticTensor(t, x, temperature);
std::tie(eps_tensile, sigma, sigma_tensile, sigma_compressive,
D, strain_energy_tensile, elastic_energy) =
MaterialLib::Solids::Phasefield::
calculateOrthoVolDevDegradedStress<DisplacementDim>(
degradation, eps, C_ortho);
break;
}
case EnergySplitModel::OrthoMasonry:
{
double temperature = 0.;
MathLib::KelvinVector::KelvinMatrixType<DisplacementDim>
C_ortho = static_cast<
MaterialLib::Solids::LinearElasticOrthotropic<
DisplacementDim> const&>(solid_material)
.getElasticTensor(t, x, temperature);
std::tie(eps_tensile, sigma, sigma_tensile, D,
strain_energy_tensile, elastic_energy) =
MaterialLib::Solids::Phasefield::
calculateOrthoMasonryDegradedStress<DisplacementDim>(
degradation, eps, C_ortho);
break;
}
}
history_variable =
......@@ -188,13 +221,6 @@ public:
_process_data.solid_materials,
_process_data.material_ids,
e.getID());
if (!dynamic_cast<MaterialLib::Solids::LinearElasticIsotropic<
DisplacementDim> const*>(&solid_material))
{
OGS_FATAL(
"For phasefield process only linear elastic solid material "
"support is implemented.");
}
auto const shape_matrices =
NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
......@@ -215,11 +241,11 @@ public:
static const int kelvin_vector_size =
MathLib::KelvinVector::kelvin_vector_dimensions(
DisplacementDim);
ip_data.eps_tensile.setZero(kelvin_vector_size);
ip_data.eps.setZero(kelvin_vector_size);
ip_data.eps_prev.resize(kelvin_vector_size);
ip_data.C_tensile.setZero(kelvin_vector_size, kelvin_vector_size);
ip_data.C_compressive.setZero(kelvin_vector_size,
kelvin_vector_size);
ip_data.D.setZero(kelvin_vector_size, kelvin_vector_size);
ip_data.sigma_tensile.setZero(kelvin_vector_size);
ip_data.sigma_compressive.setZero(kelvin_vector_size);
ip_data.sigma.setZero(kelvin_vector_size);
......@@ -313,6 +339,26 @@ private:
_ip_data, &IpData::sigma, cache);
}
std::vector<double> const& getIntPtSigmaTensile(
const double /*t*/,
std::vector<GlobalVector*> const& /*x*/,
std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
std::vector<double>& cache) const override
{
return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
_ip_data, &IpData::sigma_tensile, cache);
}
std::vector<double> const& getIntPtSigmaCompressive(
const double /*t*/,
std::vector<GlobalVector*> const& /*x*/,
std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
std::vector<double>& cache) const override
{
return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
_ip_data, &IpData::sigma_compressive, cache);
}
std::vector<double> const& getIntPtEpsilon(
const double /*t*/,
std::vector<GlobalVector*> const& /*x*/,
......@@ -323,6 +369,16 @@ private:
_ip_data, &IpData::eps, cache);
}
std::vector<double> const& getIntPtEpsilonTensile(
const double /*t*/,
std::vector<GlobalVector*> const& /*x*/,
std::vector<NumLib::LocalToGlobalIndexMap const*> const& /*dof_table*/,
std::vector<double>& cache) const override
{
return ProcessLib::getIntegrationPointKelvinVectorData<DisplacementDim>(
_ip_data, &IpData::eps_tensile, cache);
}
void assembleWithJacobianPhaseFieldEquations(
double const t, double const dt, Eigen::VectorXd const& local_x,
std::vector<double>& local_b_data, std::vector<double>& local_Jac_data);
......
......@@ -137,6 +137,27 @@ void PhaseFieldProcess<DisplacementDim>::initializeConcreteProcess(
getExtrapolator(), _local_assemblers,
&LocalAssemblerInterface::getIntPtEpsilon));
_secondary_variables.addSecondaryVariable(
"sigma_tensile",
makeExtrapolator(MathLib::KelvinVector::KelvinVectorType<
DisplacementDim>::RowsAtCompileTime,
getExtrapolator(), _local_assemblers,
&LocalAssemblerInterface::getIntPtSigmaTensile));
_secondary_variables.addSecondaryVariable(
"sigma_compressive",
makeExtrapolator(MathLib::KelvinVector::KelvinVectorType<
DisplacementDim>::RowsAtCompileTime,
getExtrapolator(), _local_assemblers,
&LocalAssemblerInterface::getIntPtSigmaCompressive));
_secondary_variables.addSecondaryVariable(
"eps_tensile",
makeExtrapolator(MathLib::KelvinVector::KelvinVectorType<
DisplacementDim>::RowsAtCompileTime,
getExtrapolator(), _local_assemblers,
&LocalAssemblerInterface::getIntPtEpsilonTensile));
// Initialize local assemblers after all variables have been set.
GlobalExecutor::executeMemberOnDereferenced(
&LocalAssemblerInterface::initialize, _local_assemblers,
......
......@@ -49,7 +49,9 @@ enum class EnergySplitModel
{
Isotropic,
VolDev,
EffectiveStress
EffectiveStress,
OrthoVolDev,
OrthoMasonry
};
class DegradationDerivative
......
......@@ -119,3 +119,18 @@ AddTest(
expected_2D_PropagatingCrack_AT1_h0p01_ts_2_t_0.020000.vtu 2D_PropagatingCrack_AT1_h0p01_ts_2_t_0.020000.vtu displacement displacement 1e-5 0
expected_2D_PropagatingCrack_AT1_h0p01_ts_2_t_0.020000.vtu 2D_PropagatingCrack_AT1_h0p01_ts_2_t_0.020000.vtu phasefield phasefield 1e-6 0
)
AddTest(
NAME PhaseField_3D_beam_tens_AT2_vd_ortho
PATH PhaseField/beam/voldev-ortho
EXECUTABLE ogs
EXECUTABLE_ARGS AT2_vd_tensile_VZ.prj
WRAPPER mpirun
WRAPPER_ARGS -np 1
TESTER vtkdiff
REQUIREMENTS OGS_USE_MPI
RUNTIME 120
DIFF_DATA
expected_AT2_vd_tensile_ts_10_t_1.000000.vtu AT2_vd_tensile_ts_10_t_1.000000.vtu displacement displacement 1e-5 0
expected_AT2_vd_tensile_ts_10_t_1.000000.vtu AT2_vd_tensile_ts_10_t_1.000000.vtu phasefield phasefield 1e-6 0
)
\ No newline at end of file
......@@ -10,10 +10,9 @@
*/
#pragma once
#include "ThermoMechanicalPhaseFieldFEM.h"
#include "NumLib/DOF/DOFTableUtil.h"
#include "ProcessLib/CoupledSolutionsForStaggeredScheme.h"
#include "ThermoMechanicalPhaseFieldFEM.h"
namespace ProcessLib
{
......@@ -93,9 +92,7 @@ void ThermoMechanicalPhaseFieldLocalAssembler<ShapeFunction, DisplacementDim>::
auto& eps = _ip_data[ip].eps;
auto const& C_tensile = _ip_data[ip].C_tensile;
auto const& C_compressive = _ip_data[ip].C_compressive;
auto const& D = _ip_data[ip].D;
auto const& sigma = _ip_data[ip].sigma;
auto const k = _process_data.residual_stiffness(t, x_position)[0];
......@@ -113,7 +110,6 @@ void ThermoMechanicalPhaseFieldLocalAssembler<ShapeFunction, DisplacementDim>::
double const d_ip = N.dot(d);
double const degradation = d_ip * d_ip * (1 - k) + k;
auto const C_eff = degradation * C_tensile + C_compressive;
eps.noalias() = B * u;
_ip_data[ip].updateConstitutiveRelation(t, x_position, dt, u, alpha,
delta_T, degradation);
......@@ -131,7 +127,7 @@ void ThermoMechanicalPhaseFieldLocalAssembler<ShapeFunction, DisplacementDim>::
.noalias() = N;
}
local_Jac.noalias() += B.transpose() * C_eff * B * w;
local_Jac.noalias() += B.transpose() * D * B * w;
local_rhs.noalias() -=
(B.transpose() * sigma - N_u.transpose() * rho_s * b) * w;
......
......@@ -60,7 +60,7 @@ struct IntegrationPointData final
DisplacementDim>::MaterialStateVariables>
material_state_variables;
typename BMatricesType::KelvinMatrixType C_tensile, C_compressive;
typename BMatricesType::KelvinMatrixType D;
double integration_weight;
void pushBackState()
......@@ -91,8 +91,8 @@ struct IntegrationPointData final
auto const bulk_modulus = linear_elastic_mp.bulk_modulus(t, x);
auto const mu = linear_elastic_mp.mu(t, x);
std::tie(sigma, sigma_tensile, C_tensile, C_compressive,
strain_energy_tensile, elastic_energy) =
std::tie(sigma, sigma_tensile, D, strain_energy_tensile,
elastic_energy) =
MaterialLib::Solids::Phasefield::calculateVolDevDegradedStress<
DisplacementDim>(degradation, bulk_modulus, mu, eps_m);
}
......@@ -206,9 +206,8 @@ public:
ip_data.eps.setZero(kelvin_vector_size);
ip_data.eps_prev.resize(kelvin_vector_size);
ip_data.eps_m.setZero(kelvin_vector_size);
ip_data.C_tensile.setZero(kelvin_vector_size, kelvin_vector_size);
ip_data.C_compressive.setZero(kelvin_vector_size,
kelvin_vector_size);
ip_data.D.setZero(kelvin_vector_size, kelvin_vector_size);
ip_data.sigma_tensile.setZero(kelvin_vector_size);
ip_data.sigma_compressive.setZero(kelvin_vector_size);
ip_data.heatflux.setZero(DisplacementDim);
......
<?xml version="1.0" encoding="ISO-8859-1"?>
<OpenGeoSysProject>
<meshes>
<mesh>bar.vtu</mesh>
<mesh>bar_left.vtu</mesh>
<mesh>bar_right.vtu</mesh>
</meshes>
<processes>
<process>
<name>PhaseField</name>
<type>PHASE_FIELD</type>
<coupling_scheme>staggered</coupling_scheme>
<integration_order>2</integration_order>
<phasefield_model>AT2</phasefield_model>
<energy_split_model>OrthoVolDev</energy_split_model>
<irreversible_threshold>0.05</irreversible_threshold>
<constitutive_relation>
<type>LinearElasticOrthotropic</type>
<youngs_moduli>E</youngs_moduli>
<poissons_ratios>nu</poissons_ratios>
<shear_moduli>G</shear_moduli>
</constitutive_relation>
<phasefield_parameters>
<residual_stiffness>k</residual_stiffness>
<crack_resistance>gc</crack_resistance>
<crack_length_scale>ls</crack_length_scale>
</phasefield_parameters>
<solid_density>rho_sr</solid_density>
<process_variables>
<phasefield>phasefield</phasefield>
<displacement>displacement</displacement>
</process_variables>
<secondary_variables>
<secondary_variable type="static" internal_name="sigma" output_name="sigma"/>
<secondary_variable type="static" internal_name="sigma_tensile" output_name="sigma_tensile"/>
<secondary_variable type="static" internal_name="epsilon" output_name="epsilon"/>
<secondary_variable type="static" internal_name="eps_tensile" output_name="eps_tensile"/>
</secondary_variables>
<specific_body_force>0 0 0</specific_body_force>
</process>
</processes>
<time_loop>
<global_process_coupling>
<max_iter> 1000 </max_iter>
<convergence_criteria>
<!-- convergence criterion for the first process -->
<convergence_criterion>
<type>DeltaX</type>
<norm_type>INFINITY_N</norm_type>
<abstol>1.e-1</abstol>
<reltol>1.e-1</reltol>
</convergence_criterion>
<!-- convergence criterion for the second process -->
<convergence_criterion>
<type>DeltaX</type>
<norm_type>INFINITY_N</norm_type>
<abstol>1.e-4</abstol>
<reltol>1.e-10</reltol>
</convergence_criterion>
</convergence_criteria>
</global_process_coupling>
<processes>
<process ref="PhaseField">
<nonlinear_solver>basic_newton_u</nonlinear_solver>
<convergence_criterion>
<type>DeltaX</type>
<norm_type>NORM2</norm_type>
<reltol>1.e-14</reltol>
</convergence_criterion>
<time_discretization>
<type>BackwardEuler</type>
</time_discretization>
<time_stepping>
<type>FixedTimeStepping</type>
<t_initial>0</t_initial>
<t_end>1.0</t_end>
<timesteps>
<pair>
<repeat>1</repeat>
<delta_t>1.e-1</delta_t>
</pair>
</timesteps>
</time_stepping>
</process>
<process ref="PhaseField">
<nonlinear_solver>petsc_snes</nonlinear_solver>
<convergence_criterion>
<type>DeltaX</type>
<norm_type>NORM2</norm_type>
<reltol>1.e-14</reltol>
</convergence_criterion>
<time_discretization>
<type>BackwardEuler</type>
</time_discretization>
<time_stepping>
<type>FixedTimeStepping</type>
<t_initial>0</t_initial>
<t_end>1.0</t_end>
<timesteps>
<pair>
<repeat>1</repeat>
<delta_t>1.e-1</delta_t>
</pair>
</timesteps>
</time_stepping>
</process>
</processes>
<output>
<variables>
<variable>displacement</variable>
<variable>phasefield</variable>
<variable>sigma</variable>
<variable>sigma_tensile</variable>
<variable>epsilon</variable>
<variable>epsilon_tensile</variable>
</variables>
<type>VTK</type>
<prefix>AT2_vd_tensile</prefix>
<timesteps>
<pair>
<repeat>10</repeat>
<each_steps>1</each_steps>
</pair>
<pair>
<repeat>10</repeat>
<each_steps>1</each_steps>
</pair>
</timesteps>
</output>
</time_loop>
<local_coordinate_system>
<basis_vector_0>e0</basis_vector_0>
<basis_vector_1>e1</basis_vector_1>
<basis_vector_2>e2</basis_vector_2>
</local_coordinate_system>
<parameters>
<!-- Mechanics -->
<parameter>
<name>E</name>
<type>Constant</type>
<values>1.0 1.0 1.0</values>
</parameter>
<parameter>
<name>nu</name>
<type>Constant</type>
<values>0.15 0.15 0.15</values>
</parameter>
<parameter>
<name>G</name>
<type>Constant</type>
<values>0.4347826086956522 0.4347826086956522 0.4347826086956522</values>
</parameter>
<parameter>
<name>gc</name>
<type>Constant</type>
<value>1.0</value>
</parameter>
<parameter>
<name>e0</name>
<type>Constant</type>
<values>0.7071067811865475 0.7071067811865475 0.0</values>
</parameter>
<parameter>
<name>e1</name>
<type>Constant</type>
<values>-0.7071067811865475 0.7071067811865475 0.0</values>
</parameter>
<parameter>
<name>e2</name>
<type>Constant</type>
<values>0.0 0.0 1.0</values>
</parameter>
<parameter>
<name>k</name>
<type>Constant</type>
<value>1e-8</value>
</parameter>
<parameter>
<name>ls</name>
<type>Constant</type>
<value>0.02</value>
</parameter>
<parameter>
<name>rho_sr</name>
<type>Constant</type>
<value>0.0</value>
</parameter>
<parameter>
<name>displacement0</name>
<type>Constant</type>
<values>0 0 0</values>
</parameter>
<parameter>
<name>phasefield_ic</name>
<type>Constant</type>
<value>1</value>
</parameter>
<parameter>
<name>phasefield_bc</name>
<type>Constant</type>
<value>1</value>
</parameter>
<parameter>
<name>dirichlet0</name>
<type>Constant</type>
<value>0</value>
</parameter>
<parameter>
<name>Dirichlet_spatial</name>
<type>Constant</type>
<value>1</value>
</parameter>
<parameter>
<name>dirichlet_right_time</name>
<type>CurveScaled</type>
<curve>dirichlet_time</curve>
<parameter>dirichlet_right</parameter>
</parameter>
<parameter>
<name>dirichlet_right</name>
<type>Constant</type>
<value>5.0</value>
</parameter>
</parameters>
<curves>
<curve>
<name>dirichlet_time</name>
<coords>0 1.0</coords>
<values>0 0.9</values>
</curve>
</curves>
<process_variables>
<process_variable>
<name>phasefield</name>
<components>1</components>
<order>1</order>
<initial_condition>phasefield_ic</initial_condition>
<boundary_conditions>
<boundary_condition>
<geometrical_set>bar</geometrical_set>
<geometry>left</geometry>
<type>Dirichlet</type>
<component>0</component>
<parameter>phasefield_bc</parameter>
</boundary_condition>
<boundary_condition>
<geometrical_set>bar</geometrical_set>
<geometry>right</geometry>
<type>Dirichlet</type>
<component>0</component>
<parameter>phasefield_bc</parameter>
</boundary_condition>
<boundary_condition>
<geometrical_set>bar</geometrical_set>
<geometry>right</geometry>
<type>PhaseFieldIrreversibleDamageOracleBoundaryCondition</type>
<component>0</component>
</boundary_condition>
</boundary_conditions>
</process_variable>
<process_variable>
<name>displacement</name>
<components>3</components>
<order>1</order>
<initial_condition>displacement0</initial_condition>
<boundary_conditions>
<boundary_condition>
<geometrical_set>bar</geometrical_set>
<geometry>left</geometry>
<type>Dirichlet</type>
<component>0</component>
<parameter>dirichlet0</parameter>
</boundary_condition>
<boundary_condition>
<geometrical_set>bar</geometrical_set>
<geometry>right</geometry>
<type>Dirichlet</type>
<component>0</component>
<parameter>dirichlet_right_time</parameter>
</boundary_condition>
</boundary_conditions>
</process_variable>
</process_variables>
<nonlinear_solvers>
<nonlinear_solver>
<name>petsc_snes</name>
<type>PETScSNES</type>
<max_iter>50</max_iter>
<linear_solver>linear_solver_d</linear_solver>
</nonlinear_solver>
<nonlinear_solver>
<name>basic_newton_u</name>
<type>Newton</type>
<max_iter>200</max_iter>
<linear_solver>linear_solver_u</linear_solver>
</nonlinear_solver>
</nonlinear_solvers>
<linear_solvers>
<linear_solver>
<name>linear_solver_d</name>
<eigen>
<solver_type>BiCGSTAB</solver_type>
<precon_type>ILUT</precon_type>
<max_iteration_step>10000</max_iteration_step>
<error_tolerance>1e-16</error_tolerance>
</eigen>
<petsc>
<parameters>-ksp_type cg -pc_type bjacobi -ksp_atol 1e-10 -ksp_rtol 1e-10 -snes_type vinewtonrsls -snes_linesearch_type l2 -snes_atol 1.e-8 -snes_rtol 1.e-8 -snes_max_it 1000 -snes_monitor </parameters>
</petsc>
</linear_solver>
<linear_solver>
<name>linear_solver_u</name>
<petsc>
<parameters>-ksp_type cg -pc_type bjacobi -ksp_atol 1e-14 -ksp_rtol 1e-14 </parameters>
</petsc>
</linear_solver>
</linear_solvers>
</OpenGeoSysProject>
<?xml version="1.0" encoding="ISO-8859-1"?>
<?xml-stylesheet type="text/xsl" href="OpenGeoSysGLI.xsl"?>
<OpenGeoSysGLI xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:noNamespaceSchemaLocation="http://www.opengeosys.org/images/xsd/OpenGeoSysGLI.xsd" xmlns:ogs="http://www.opengeosys.org">
<name>bar</name>
<points>
<point id="0" x="0" y="0" z="0" name="p_0"/>
<point id="1" x="1.0" y="0" z="0" name="p_1"/>
<point id="2" x="1.0" y="0.05" z="0" name="p_2"/>
<point id="3" x="0" y="0.05" z="0" name="p_3"/>
<point id="4" x="0" y="0" z="0.05" name="p_4"/>
<point id="5" x="1.0" y="0" z="0.05" name="p_5"/>
<point id="6" x="1.0" y="0.05" z="0.05" name="p_6"/>
<point id="7" x="0" y="0.05" z="0.05" name="p_7"/>
</points>
<polylines>
<polyline id="0" name="left_bot">
<pnt>0</pnt>
<pnt>4</pnt>
</polyline>
<polyline id="1" name="left_front">
<pnt>0</pnt>
<pnt>3</pnt>
</polyline>
</polylines>
<surfaces>
<surface id="0" name="front">
<element p1="0" p2="1" p3="3"/>
<element p1="1" p2="3" p3="2"/>
</surface>
<surface id="1" name="back">
<element p1="4" p2="5" p3="7"/>
<element p1="5" p2="7" p3="6"/>
</surface>
<surface id="2" name="left">
<element p1="0" p2="4" p3="3"/>
<element p1="4" p2="3" p3="7"/>
</surface>
<surface id="3" name="right">
<element p1="1" p2="5" p3="2"/>
<element p1="5" p2="2" p3="6"/>
</surface>
<surface id="4" name="top">
<element p1="3" p2="2" p3="7"/>
<element p1="2" p2="7" p3="6"/>
</surface>
<surface id="5" name="bottom">
<element p1="0" p2="1" p3="4"/>
<element p1="1" p2="4" p3="5"/>
</surface>
</surfaces>
</OpenGeoSysGLI>
File added
<?xml version="1.0"?>
<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64">
<UnstructuredGrid>
<Piece NumberOfPoints="36" NumberOfCells="25" >
<PointData>
<DataArray type="UInt64" Name="bulk_node_ids" format="appended" RangeMin="0" RangeMax="71" offset="0" />
</PointData>
<CellData>
<DataArray type="UInt64" Name="bulk_element_ids" format="appended" RangeMin="0" RangeMax="2400" offset="396" />
</CellData>
<Points>
<DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0" RangeMax="0.070710678119" offset="676" />
</Points>
<Cells>
<DataArray type="Int64" Name="connectivity" format="appended" RangeMin="" RangeMax="" offset="1840" />
<DataArray type="Int64" Name="offsets" format="appended" RangeMin="" RangeMax="" offset="2920" />
<DataArray type="UInt8" Name="types" format="appended" RangeMin="" RangeMax="" offset="3200" />
</Cells>
</Piece>
</UnstructuredGrid>
<AppendedData encoding="base64">
_IAEAAAAAAAAXAAAAAAAAABYAAAAAAAAAFQAAAAAAAAAUAAAAAAAAADkAAAAAAAAAAAAAAAAAAAAIAAAAAAAAADgAAAAAAAAAQQAAAAAAAAATAAAAAAAAAD8AAAAAAAAACgAAAAAAAAA+AAAAAAAAAD0AAAAAAAAACQAAAAAAAAASAAAAAAAAAEIAAAAAAAAAAwAAAAAAAABAAAAAAAAAAEcAAAAAAAAADwAAAAAAAAA6AAAAAAAAAA4AAAAAAAAAOwAAAAAAAAANAAAAAAAAAEMAAAAAAAAAAQAAAAAAAAAQAAAAAAAAAEUAAAAAAAAAAgAAAAAAAAARAAAAAAAAAAsAAAAAAAAARAAAAAAAAAA8AAAAAAAAAEYAAAAAAAAADAAAAAAAAAA=yAAAAAAAAAAAAAAAAAAAAGQAAAAAAAAAyAAAAAAAAAAsAQAAAAAAAJABAAAAAAAA9AEAAAAAAABYAgAAAAAAALwCAAAAAAAAIAMAAAAAAACEAwAAAAAAAOgDAAAAAAAATAQAAAAAAACwBAAAAAAAABQFAAAAAAAAeAUAAAAAAADcBQAAAAAAAEAGAAAAAAAApAYAAAAAAAAIBwAAAAAAAGwHAAAAAAAA0AcAAAAAAAA0CAAAAAAAAJgIAAAAAAAA/AgAAAAAAABgCQAAAAAAAA==YAMAAAAAAAAAAAAAAAAAAHsUrkfheqQ/AAAAAAAAAAAAAAAAAAAAALgehetRuJ4/AAAAAAAAAAAAAAAAAAAAAHsUrkfhepQ/AAAAAAAAAAAAAAAAAAAAAHsUrkfheoQ/AAAAAAAAAAAAAAAAAAAAAHsUrkfhepQ/exSuR+F6hD8AAAAAAAAAAAAAAAAAAAAAmpmZmZmZqT8AAAAAAAAAAAAAAAAAAAAAexSuR+F6hD8AAAAAAAAAAHsUrkfheoQ/fRSuR+F6hD8AAAAAAAAAAHwUrkfhepQ/uB6F61G4nj8AAAAAAAAAAJqZmZmZmak/exSuR+F6pD8AAAAAAAAAAHsUrkfheqQ/fBSuR+F6lD8AAAAAAAAAAAAAAAAAAAAAuB6F61G4nj8AAAAAAAAAALgehetRuJ4/fBSuR+F6lD8AAAAAAAAAAHsUrkfhepQ/exSuR+F6lD8AAAAAAAAAAAAAAAAAAAAAexSuR+F6lD8AAAAAAAAAAJqZmZmZmak/uB6F61G4nj8AAAAAAAAAALgehetRuJ4/uB6F61G4nj8AAAAAAAAAAJqZmZmZmak/AAAAAAAAAAAAAAAAAAAAAHwUrkfheoQ/uB6F61G4nj8AAAAAAAAAAH4UrkfheqQ/fBSuR+F6pD8AAAAAAAAAAHsUrkfheqQ/mpmZmZmZqT8AAAAAAAAAALgehetRuJ4/fBSuR+F6hD8AAAAAAAAAALgehetRuJ4/mpmZmZmZqT8AAAAAAAAAAHsUrkfheqQ/fBSuR+F6hD8AAAAAAAAAAHsUrkfhepQ/mpmZmZmZqT8AAAAAAAAAAHwUrkfheqQ/uB6F61G4nj8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAJqZmZmZmak/exSuR+F6hD8AAAAAAAAAAH4UrkfhepQ/exSuR+F6pD8AAAAAAAAAAJqZmZmZmak/mpmZmZmZqT8AAAAAAAAAAJqZmZmZmak/exSuR+F6lD8AAAAAAAAAAAAAAAAAAAAAexSuR+F6pD8AAAAAAAAAAH4UrkfheoQ/fRSuR+F6pD8AAAAAAAAAAHsUrkfheoQ/fRSuR+F6lD8AAAAAAAAAALgehetRuJ4/fBSuR+F6pD8AAAAAAAAAAHsUrkfheoQ/mpmZmZmZqT8=IAMAAAAAAAAFAAAAAAAAACMAAAAAAAAAIAAAAAAAAAAfAAAAAAAAACMAAAAAAAAAGAAAAAAAAAAcAAAAAAAAACAAAAAAAAAAGAAAAAAAAAAWAAAAAAAAACIAAAAAAAAAHAAAAAAAAAAWAAAAAAAAABQAAAAAAAAAEwAAAAAAAAAiAAAAAAAAABQAAAAAAAAAHQAAAAAAAAAJAAAAAAAAABMAAAAAAAAAHwAAAAAAAAAgAAAAAAAAABIAAAAAAAAACwAAAAAAAAAgAAAAAAAAABwAAAAAAAAACAAAAAAAAAASAAAAAAAAABwAAAAAAAAAIgAAAAAAAAAQAAAAAAAAAAgAAAAAAAAAIgAAAAAAAAATAAAAAAAAABkAAAAAAAAAEAAAAAAAAAATAAAAAAAAAAkAAAAAAAAADwAAAAAAAAAZAAAAAAAAAAsAAAAAAAAAEgAAAAAAAAAhAAAAAAAAAA4AAAAAAAAAEgAAAAAAAAAIAAAAAAAAAA0AAAAAAAAAIQAAAAAAAAAIAAAAAAAAABAAAAAAAAAADAAAAAAAAAANAAAAAAAAABAAAAAAAAAAGQAAAAAAAAAKAAAAAAAAAAwAAAAAAAAAGQAAAAAAAAAPAAAAAAAAAB4AAAAAAAAACgAAAAAAAAAOAAAAAAAAACEAAAAAAAAABwAAAAAAAAAGAAAAAAAAACEAAAAAAAAADQAAAAAAAAAEAAAAAAAAAAcAAAAAAAAADQAAAAAAAAAMAAAAAAAAABUAAAAAAAAABAAAAAAAAAAMAAAAAAAAAAoAAAAAAAAAFwAAAAAAAAAVAAAAAAAAAAoAAAAAAAAAHgAAAAAAAAAbAAAAAAAAABcAAAAAAAAABgAAAAAAAAAHAAAAAAAAAAMAAAAAAAAAGgAAAAAAAAAHAAAAAAAAAAQAAAAAAAAAAgAAAAAAAAADAAAAAAAAAAQAAAAAAAAAFQAAAAAAAAABAAAAAAAAAAIAAAAAAAAAFQAAAAAAAAAXAAAAAAAAAAAAAAAAAAAAAQAAAAAAAAAXAAAAAAAAABsAAAAAAAAAEQAAAAAAAAAAAAAAAAAAAA==yAAAAAAAAAAEAAAAAAAAAAgAAAAAAAAADAAAAAAAAAAQAAAAAAAAABQAAAAAAAAAGAAAAAAAAAAcAAAAAAAAACAAAAAAAAAAJAAAAAAAAAAoAAAAAAAAACwAAAAAAAAAMAAAAAAAAAA0AAAAAAAAADgAAAAAAAAAPAAAAAAAAABAAAAAAAAAAEQAAAAAAAAASAAAAAAAAABMAAAAAAAAAFAAAAAAAAAAVAAAAAAAAABYAAAAAAAAAFwAAAAAAAAAYAAAAAAAAABkAAAAAAAAAA==GQAAAAAAAAAJCQkJCQkJCQkJCQkJCQkJCQkJCQkJCQkJ
</AppendedData>
</VTKFile>
<?xml version="1.0"?>
<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64">
<UnstructuredGrid>
<Piece NumberOfPoints="36" NumberOfCells="25" >
<PointData>
<DataArray type="UInt64" Name="bulk_node_ids" format="appended" RangeMin="4" RangeMax="87" offset="0" />
</PointData>
<CellData>
<DataArray type="UInt64" Name="bulk_element_ids" format="appended" RangeMin="99" RangeMax="2499" offset="396" />
</CellData>
<Points>
<DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="1" RangeMax="1.0024968828" offset="676" />
</Points>
<Cells>
<DataArray type="Int64" Name="connectivity" format="appended" RangeMin="" RangeMax="" offset="1840" />
<DataArray type="Int64" Name="offsets" format="appended" RangeMin="" RangeMax="" offset="2920" />
<DataArray type="UInt8" Name="types" format="appended" RangeMin="" RangeMax="" offset="3200" />
</Cells>
</Piece>
</UnstructuredGrid>
<AppendedData encoding="base64">
_IAEAAAAAAAAHAAAAAAAAACUAAAAAAAAAJwAAAAAAAAAkAAAAAAAAAAUAAAAAAAAAJgAAAAAAAAAgAAAAAAAAAAYAAAAAAAAAVQAAAAAAAABUAAAAAAAAAFIAAAAAAAAAHwAAAAAAAABTAAAAAAAAAFEAAAAAAAAAVgAAAAAAAAAhAAAAAAAAACIAAAAAAAAAHQAAAAAAAABMAAAAAAAAAE8AAAAAAAAAGgAAAAAAAAAjAAAAAAAAAEsAAAAAAAAAHgAAAAAAAABNAAAAAAAAAEoAAAAAAAAAVwAAAAAAAABJAAAAAAAAAFAAAAAAAAAAHAAAAAAAAABIAAAAAAAAAE4AAAAAAAAAGwAAAAAAAAAYAAAAAAAAABkAAAAAAAAABAAAAAAAAAA=yAAAAAAAAABjAAAAAAAAAMcAAAAAAAAAKwEAAAAAAACPAQAAAAAAAPMBAAAAAAAAVwIAAAAAAAC7AgAAAAAAAB8DAAAAAAAAgwMAAAAAAADnAwAAAAAAAEsEAAAAAAAArwQAAAAAAAATBQAAAAAAAHcFAAAAAAAA2wUAAAAAAAA/BgAAAAAAAKMGAAAAAAAABwcAAAAAAABrBwAAAAAAAM8HAAAAAAAAMwgAAAAAAACXCAAAAAAAAPsIAAAAAAAAXwkAAAAAAADDCQAAAAAAAA==YAMAAAAAAAAAAAAAAADwP5qZmZmZmak/AAAAAAAAAAAAAAAAAADwP3sUrkfhepQ/AAAAAAAAAAAAAAAAAADwP3sUrkfheqQ/AAAAAAAAAAAAAAAAAADwP3sUrkfheoQ/AAAAAAAAAAAAAAAAAADwPwAAAAAAAAAAAAAAAAAAAAAAAAAAAADwP7gehetRuJ4/AAAAAAAAAAAAAAAAAADwP5qZmZmZmak/exSuR+F6hD8AAAAAAADwP5qZmZmZmak/mpmZmZmZqT8AAAAAAADwP34UrkfhepQ/exSuR+F6hD8AAAAAAADwP34UrkfheoQ/exSuR+F6hD8AAAAAAADwP7gehetRuJ4/fBSuR+F6lD8AAAAAAADwP3sUrkfheqQ/mpmZmZmZqT8AAAAAAADwP3wUrkfheqQ/fBSuR+F6lD8AAAAAAADwP3wUrkfhepQ/fBSuR+F6lD8AAAAAAADwP7gehetRuJ4/exSuR+F6hD8AAAAAAADwP5qZmZmZmak/exSuR+F6lD8AAAAAAADwP5qZmZmZmak/uB6F61G4nj8AAAAAAADwP3sUrkfhepQ/mpmZmZmZqT8AAAAAAADwP3wUrkfheoQ/uR6F61G4nj8AAAAAAADwP3wUrkfheqQ/uR6F61G4nj8AAAAAAADwPwAAAAAAAAAAuB6F61G4nj8AAAAAAADwP5qZmZmZmak/exSuR+F6pD8AAAAAAADwP3wUrkfheqQ/fhSuR+F6pD8AAAAAAADwP7gehetRuJ4/mpmZmZmZqT8AAAAAAADwP3wUrkfhepQ/uR6F61G4nj8AAAAAAADwP7gehetRuJ4/fBSuR+F6pD8AAAAAAADwP34UrkfheqQ/exSuR+F6hD8AAAAAAADwP3wUrkfhepQ/fRSuR+F6pD8AAAAAAADwP3wUrkfheoQ/fBSuR+F6lD8AAAAAAADwP3sUrkfheoQ/mpmZmZmZqT8AAAAAAADwP3wUrkfheoQ/exSuR+F6pD8AAAAAAADwP7cehetRuJ4/uB6F61G4nj8AAAAAAADwPwAAAAAAAAAAexSuR+F6pD8AAAAAAADwPwAAAAAAAAAAexSuR+F6hD8AAAAAAADwPwAAAAAAAAAAexSuR+F6lD8AAAAAAADwPwAAAAAAAAAAmpmZmZmZqT8=IAMAAAAAAAAjAAAAAAAAACAAAAAAAAAAHgAAAAAAAAAdAAAAAAAAAB0AAAAAAAAAHgAAAAAAAAAbAAAAAAAAABEAAAAAAAAAEQAAAAAAAAAbAAAAAAAAABkAAAAAAAAAFwAAAAAAAAAXAAAAAAAAABkAAAAAAAAAFgAAAAAAAAALAAAAAAAAAAsAAAAAAAAAFgAAAAAAAAAVAAAAAAAAAAcAAAAAAAAAIAAAAAAAAAAUAAAAAAAAABIAAAAAAAAAHgAAAAAAAAAeAAAAAAAAABIAAAAAAAAAGAAAAAAAAAAbAAAAAAAAABsAAAAAAAAAGAAAAAAAAAAfAAAAAAAAABkAAAAAAAAAGQAAAAAAAAAfAAAAAAAAABMAAAAAAAAAFgAAAAAAAAAWAAAAAAAAABMAAAAAAAAAEAAAAAAAAAAVAAAAAAAAABQAAAAAAAAAIgAAAAAAAAAcAAAAAAAAABIAAAAAAAAAEgAAAAAAAAAcAAAAAAAAAA0AAAAAAAAAGAAAAAAAAAAYAAAAAAAAAA0AAAAAAAAACgAAAAAAAAAfAAAAAAAAAB8AAAAAAAAACgAAAAAAAAAMAAAAAAAAABMAAAAAAAAAEwAAAAAAAAAMAAAAAAAAAA8AAAAAAAAAEAAAAAAAAAAiAAAAAAAAACEAAAAAAAAACQAAAAAAAAAcAAAAAAAAABwAAAAAAAAACQAAAAAAAAAIAAAAAAAAAA0AAAAAAAAADQAAAAAAAAAIAAAAAAAAAA4AAAAAAAAACgAAAAAAAAAKAAAAAAAAAA4AAAAAAAAAGgAAAAAAAAAMAAAAAAAAAAwAAAAAAAAAGgAAAAAAAAAGAAAAAAAAAA8AAAAAAAAAIQAAAAAAAAAEAAAAAAAAAAMAAAAAAAAACQAAAAAAAAAJAAAAAAAAAAMAAAAAAAAAAQAAAAAAAAAIAAAAAAAAAAgAAAAAAAAAAQAAAAAAAAAFAAAAAAAAAA4AAAAAAAAADgAAAAAAAAAFAAAAAAAAAAIAAAAAAAAAGgAAAAAAAAAaAAAAAAAAAAIAAAAAAAAAAAAAAAAAAAAGAAAAAAAAAA==yAAAAAAAAAAEAAAAAAAAAAgAAAAAAAAADAAAAAAAAAAQAAAAAAAAABQAAAAAAAAAGAAAAAAAAAAcAAAAAAAAACAAAAAAAAAAJAAAAAAAAAAoAAAAAAAAACwAAAAAAAAAMAAAAAAAAAA0AAAAAAAAADgAAAAAAAAAPAAAAAAAAABAAAAAAAAAAEQAAAAAAAAASAAAAAAAAABMAAAAAAAAAFAAAAAAAAAAVAAAAAAAAABYAAAAAAAAAFwAAAAAAAAAYAAAAAAAAABkAAAAAAAAAA==GQAAAAAAAAAJCQkJCQkJCQkJCQkJCQkJCQkJCQkJCQkJ
</AppendedData>
</VTKFile>
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment