Skip to content

Commit

Permalink
Merge branch '2nd_var_output_LIE' into 'master'
Browse files Browse the repository at this point in the history 
[LIE/HM] Improve output

See merge request ogs/ogs!4851
  • Loading branch information
@endJunction
endJunction committed Jan 19, 2024
2 parents 1a9aaaf + af985e7 commit 1e97abd
Show file tree
Hide file tree
Showing 54 changed files with 1,866 additions and 1,398 deletions.
233 changes: 56 additions & 177 deletions ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -212,89 +212,51 @@ void HydroMechanicsProcess<GlobalDim>::initializeConcreteProcess(
NumLib::IntegrationOrder{integration_order}, mesh.isAxiallySymmetric(),
_process_data);

auto mesh_prop_sigma_xx = MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh), "stress_xx",
MeshLib::MeshItemType::Cell, 1);
mesh_prop_sigma_xx->resize(mesh.getNumberOfElements());
_process_data.mesh_prop_stress_xx = mesh_prop_sigma_xx;

auto mesh_prop_sigma_yy = MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh), "stress_yy",
MeshLib::MeshItemType::Cell, 1);
mesh_prop_sigma_yy->resize(mesh.getNumberOfElements());
_process_data.mesh_prop_stress_yy = mesh_prop_sigma_yy;

auto mesh_prop_sigma_zz = MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh), "stress_zz",
MeshLib::MeshItemType::Cell, 1);
mesh_prop_sigma_zz->resize(mesh.getNumberOfElements());
_process_data.mesh_prop_stress_zz = mesh_prop_sigma_zz;

auto mesh_prop_sigma_xy = MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh), "stress_xy",
MeshLib::MeshItemType::Cell, 1);
mesh_prop_sigma_xy->resize(mesh.getNumberOfElements());
_process_data.mesh_prop_stress_xy = mesh_prop_sigma_xy;

if (GlobalDim == 3)
auto add_secondary_variable = [&](std::string const& name,
int const num_components,
auto get_ip_values_function)
{
auto mesh_prop_sigma_xz = MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh), "stress_xz",
MeshLib::MeshItemType::Cell, 1);
mesh_prop_sigma_xz->resize(mesh.getNumberOfElements());
_process_data.mesh_prop_stress_xz = mesh_prop_sigma_xz;
_secondary_variables.addSecondaryVariable(
name,
makeExtrapolator(num_components, getExtrapolator(),
_local_assemblers,
std::move(get_ip_values_function)));
};

auto mesh_prop_sigma_yz = MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh), "stress_yz",
MeshLib::MeshItemType::Cell, 1);
mesh_prop_sigma_yz->resize(mesh.getNumberOfElements());
_process_data.mesh_prop_stress_yz = mesh_prop_sigma_yz;
}
add_secondary_variable(
"sigma",
MathLib::KelvinVector::KelvinVectorType<GlobalDim>::RowsAtCompileTime,
&LocalAssemblerInterface::getIntPtSigma);

auto mesh_prop_epsilon_xx = MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh), "strain_xx",
MeshLib::MeshItemType::Cell, 1);
mesh_prop_epsilon_xx->resize(mesh.getNumberOfElements());
_process_data.mesh_prop_strain_xx = mesh_prop_epsilon_xx;

auto mesh_prop_epsilon_yy = MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh), "strain_yy",
MeshLib::MeshItemType::Cell, 1);
mesh_prop_epsilon_yy->resize(mesh.getNumberOfElements());
_process_data.mesh_prop_strain_yy = mesh_prop_epsilon_yy;

auto mesh_prop_epsilon_zz = MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh), "strain_zz",
MeshLib::MeshItemType::Cell, 1);
mesh_prop_epsilon_zz->resize(mesh.getNumberOfElements());
_process_data.mesh_prop_strain_zz = mesh_prop_epsilon_zz;

auto mesh_prop_epsilon_xy = MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh), "strain_xy",
MeshLib::MeshItemType::Cell, 1);
mesh_prop_epsilon_xy->resize(mesh.getNumberOfElements());
_process_data.mesh_prop_strain_xy = mesh_prop_epsilon_xy;

if (GlobalDim == 3)
{
auto mesh_prop_epsilon_xz = MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh), "strain_xz",
MeshLib::MeshItemType::Cell, 1);
mesh_prop_epsilon_xz->resize(mesh.getNumberOfElements());
_process_data.mesh_prop_strain_xz = mesh_prop_epsilon_xz;
add_secondary_variable(
"epsilon",
MathLib::KelvinVector::KelvinVectorType<GlobalDim>::RowsAtCompileTime,
&LocalAssemblerInterface::getIntPtEpsilon);

auto mesh_prop_epsilon_yz = MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh), "strain_yz",
MeshLib::MeshItemType::Cell, 1);
mesh_prop_epsilon_yz->resize(mesh.getNumberOfElements());
_process_data.mesh_prop_strain_yz = mesh_prop_epsilon_yz;
}
add_secondary_variable("velocity", GlobalDim,
&LocalAssemblerInterface::getIntPtDarcyVelocity);

add_secondary_variable("fracture_velocity", GlobalDim,
&LocalAssemblerInterface::getIntPtFractureVelocity);

add_secondary_variable("fracture_stress", GlobalDim,
&LocalAssemblerInterface::getIntPtFractureStress);

add_secondary_variable("fracture_aperture", 1,
&LocalAssemblerInterface::getIntPtFractureAperture);

add_secondary_variable(
"fracture_permeability", 1,
&LocalAssemblerInterface::getIntPtFracturePermeability);

_process_data.element_stresses = MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh), "sigma_avg",
MeshLib::MeshItemType::Cell,
MathLib::KelvinVector::KelvinVectorType<GlobalDim>::RowsAtCompileTime);

auto mesh_prop_velocity = MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh), "velocity",
_process_data.element_velocities = MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh), "velocity_avg",
MeshLib::MeshItemType::Cell, GlobalDim);
mesh_prop_velocity->resize(mesh.getNumberOfElements() * GlobalDim);
_process_data.mesh_prop_velocity = mesh_prop_velocity;

if (!_vec_fracture_elements.empty())
{
Expand All @@ -319,21 +281,26 @@ void HydroMechanicsProcess<GlobalDim>::initializeConcreteProcess(
(*mesh_prop_levelset)[e->getID()] = levelsets[0];
}

auto mesh_prop_w_n = MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh), "w_n",
MeshLib::MeshItemType::Cell, 1);
mesh_prop_w_n->resize(mesh.getNumberOfElements());
auto mesh_prop_w_s = MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh), "w_s",
MeshLib::MeshItemType::Cell, 1);
mesh_prop_w_s->resize(mesh.getNumberOfElements());
_process_data.mesh_prop_w_n = mesh_prop_w_n;
_process_data.mesh_prop_w_s = mesh_prop_w_s;
_process_data.element_local_jumps =
MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh), "local_jump_w_avg",
MeshLib::MeshItemType::Cell, GlobalDim);

_process_data.element_fracture_stresses =
MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh), "fracture_stress_avg",
MeshLib::MeshItemType::Cell, GlobalDim);

_process_data.element_fracture_velocities =
MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh), "fracture_velocity_avg",
MeshLib::MeshItemType::Cell, GlobalDim);

auto mesh_prop_b = MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh), "aperture",
const_cast<MeshLib::Mesh&>(mesh), "fracture_aperture_avg",
MeshLib::MeshItemType::Cell, 1);
mesh_prop_b->resize(mesh.getNumberOfElements());

auto const mesh_prop_matid = materialIDs(mesh);
if (!mesh_prop_matid)
{
Expand All @@ -360,27 +327,11 @@ void HydroMechanicsProcess<GlobalDim>::initializeConcreteProcess(
_process_data.mesh_prop_b = mesh_prop_b;

auto mesh_prop_k_f = MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh), "k_f",
const_cast<MeshLib::Mesh&>(mesh), "fracture_permeability_avg",
MeshLib::MeshItemType::Cell, 1);
mesh_prop_k_f->resize(mesh.getNumberOfElements());
_process_data.mesh_prop_k_f = mesh_prop_k_f;

auto mesh_prop_fracture_stress_shear =
MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh), "f_stress_s",
MeshLib::MeshItemType::Cell, 1);
mesh_prop_fracture_stress_shear->resize(mesh.getNumberOfElements());
_process_data.mesh_prop_fracture_stress_shear =
mesh_prop_fracture_stress_shear;

auto mesh_prop_fracture_stress_normal =
MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh), "f_stress_n",
MeshLib::MeshItemType::Cell, 1);
mesh_prop_fracture_stress_normal->resize(mesh.getNumberOfElements());
_process_data.mesh_prop_fracture_stress_normal =
mesh_prop_fracture_stress_normal;

auto mesh_prop_fracture_shear_failure =
MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh), "f_shear_failure",
Expand All @@ -389,36 +340,6 @@ void HydroMechanicsProcess<GlobalDim>::initializeConcreteProcess(
_process_data.mesh_prop_fracture_shear_failure =
mesh_prop_fracture_shear_failure;

auto mesh_prop_nodal_w = MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh), "nodal_w",
MeshLib::MeshItemType::Node, GlobalDim);
mesh_prop_nodal_w->resize(mesh.getNumberOfNodes() * GlobalDim);
_process_data.mesh_prop_nodal_w = mesh_prop_nodal_w;

auto mesh_prop_nodal_b = MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh), "nodal_aperture",
MeshLib::MeshItemType::Node, 1);
mesh_prop_nodal_b->resize(mesh.getNumberOfNodes());
_process_data.mesh_prop_nodal_b = mesh_prop_nodal_b;

if (GlobalDim == 3)
{
auto mesh_prop_w_s2 = MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh), "w_s2",
MeshLib::MeshItemType::Cell, 1);
mesh_prop_w_s2->resize(mesh.getNumberOfElements());
_process_data.mesh_prop_w_s2 = mesh_prop_w_s2;

auto mesh_prop_fracture_stress_shear2 =
MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh), "f_stress_s2",
MeshLib::MeshItemType::Cell, 1);
mesh_prop_fracture_stress_shear2->resize(
mesh.getNumberOfElements());
_process_data.mesh_prop_fracture_stress_shear2 =
mesh_prop_fracture_stress_shear2;
}

auto mesh_prop_nodal_p = MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh), "pressure_interpolated",
MeshLib::MeshItemType::Node, 1);
Expand Down Expand Up @@ -521,48 +442,6 @@ void HydroMechanicsProcess<GlobalDim>::postTimestepConcreteProcess(
(*x[process_id])[global_index];
}
}

// compute nodal w and aperture
auto const& R = _process_data.fracture_property->R;
auto const& b0 = _process_data.fracture_property->aperture0;
MeshLib::PropertyVector<double>& vec_w = *_process_data.mesh_prop_nodal_w;
MeshLib::PropertyVector<double>& vec_b = *_process_data.mesh_prop_nodal_b;

auto compute_nodal_aperture =
[&](std::size_t const node_id, double const w_n)
{
// skip aperture computation for element-wise defined b0 because there
// are jumps on the nodes between the element's values.
if (dynamic_cast<ParameterLib::MeshElementParameter<double> const*>(
&b0))
{
return std::numeric_limits<double>::quiet_NaN();
}

ParameterLib::SpatialPosition x;
x.setNodeID(node_id);
return w_n + b0(/*time independent*/ 0, x)[0];
};

Eigen::VectorXd g(GlobalDim);
Eigen::VectorXd w(GlobalDim);
for (MeshLib::Node const* node : _vec_fracture_nodes)
{
auto const node_id = node->getID();
g.setZero();
for (int k = 0; k < GlobalDim; k++)
{
g[k] = mesh_prop_g[node_id * GlobalDim + k];
}

w.noalias() = R * g;
for (int k = 0; k < GlobalDim; k++)
{
vec_w[node_id * GlobalDim + k] = w[k];
}

vec_b[node_id] = compute_nodal_aperture(node_id, w[GlobalDim - 1]);
}
}

template <int GlobalDim>
Expand Down
27 changes: 6 additions & 21 deletions ProcessLib/LIE/HydroMechanics/HydroMechanicsProcessData.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -110,30 +110,15 @@ struct HydroMechanicsProcessData
ParameterLib::Parameter<double> const* p0 = nullptr;

// mesh properties for output
MeshLib::PropertyVector<double>* mesh_prop_stress_xx = nullptr;
MeshLib::PropertyVector<double>* mesh_prop_stress_yy = nullptr;
MeshLib::PropertyVector<double>* mesh_prop_stress_zz = nullptr;
MeshLib::PropertyVector<double>* mesh_prop_stress_xy = nullptr;
MeshLib::PropertyVector<double>* mesh_prop_stress_yz = nullptr;
MeshLib::PropertyVector<double>* mesh_prop_stress_xz = nullptr;
MeshLib::PropertyVector<double>* mesh_prop_strain_xx = nullptr;
MeshLib::PropertyVector<double>* mesh_prop_strain_yy = nullptr;
MeshLib::PropertyVector<double>* mesh_prop_strain_zz = nullptr;
MeshLib::PropertyVector<double>* mesh_prop_strain_xy = nullptr;
MeshLib::PropertyVector<double>* mesh_prop_strain_yz = nullptr;
MeshLib::PropertyVector<double>* mesh_prop_strain_xz = nullptr;
MeshLib::PropertyVector<double>* mesh_prop_velocity = nullptr;
MeshLib::PropertyVector<double>* element_stresses = nullptr;
MeshLib::PropertyVector<double>* element_velocities = nullptr;
MeshLib::PropertyVector<double>* element_local_jumps = nullptr;
MeshLib::PropertyVector<double>* element_fracture_stresses = nullptr;
MeshLib::PropertyVector<double>* element_fracture_velocities = nullptr;

MeshLib::PropertyVector<double>* mesh_prop_b = nullptr;
MeshLib::PropertyVector<double>* mesh_prop_k_f = nullptr;
MeshLib::PropertyVector<double>* mesh_prop_w_n = nullptr;
MeshLib::PropertyVector<double>* mesh_prop_w_s = nullptr;
MeshLib::PropertyVector<double>* mesh_prop_w_s2 = nullptr;
MeshLib::PropertyVector<double>* mesh_prop_fracture_stress_shear = nullptr;
MeshLib::PropertyVector<double>* mesh_prop_fracture_stress_shear2 = nullptr;
MeshLib::PropertyVector<double>* mesh_prop_fracture_stress_normal = nullptr;
MeshLib::PropertyVector<double>* mesh_prop_fracture_shear_failure = nullptr;
MeshLib::PropertyVector<double>* mesh_prop_nodal_w = nullptr;
MeshLib::PropertyVector<double>* mesh_prop_nodal_b = nullptr;
MeshLib::PropertyVector<double>* mesh_prop_nodal_p = nullptr;

MeshLib::PropertyVector<double>* mesh_prop_nodal_forces = nullptr;
Expand Down
Loading

0 comments on commit 1e97abd

Please sign in to comment.