diff --git a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp
index 7315253f1a9ba6c10ca509e56b38f6ff7695bd77..814aa30eab88bfc37f40aae7c3a18e0f011a5469 100644
--- a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp
+++ b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp
@@ -207,66 +207,53 @@ void HydroMechanicsProcess<GlobalDim>::initializeConcreteProcess(
_local_assemblers, mesh.isAxiallySymmetric(), integration_order,
_process_data);
- auto mesh_prop_sigma_xx = const_cast<MeshLib::Mesh&>(mesh)
- .getProperties()
- .template createNewPropertyVector<double>(
- "stress_xx", MeshLib::MeshItemType::Cell);
+ 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 = const_cast<MeshLib::Mesh&>(mesh)
- .getProperties()
- .template createNewPropertyVector<double>(
- "stress_yy", MeshLib::MeshItemType::Cell);
+ 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_xy = const_cast<MeshLib::Mesh&>(mesh)
- .getProperties()
- .template createNewPropertyVector<double>(
- "stress_xy", MeshLib::MeshItemType::Cell);
+ 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;
- auto mesh_prop_epsilon_xx =
- const_cast<MeshLib::Mesh&>(mesh)
- .getProperties()
- .template createNewPropertyVector<double>(
- "strain_xx", MeshLib::MeshItemType::Cell);
+ 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 =
- const_cast<MeshLib::Mesh&>(mesh)
- .getProperties()
- .template createNewPropertyVector<double>(
- "strain_yy", MeshLib::MeshItemType::Cell);
+ 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_xy =
- const_cast<MeshLib::Mesh&>(mesh)
- .getProperties()
- .template createNewPropertyVector<double>(
- "strain_xy", MeshLib::MeshItemType::Cell);
+ 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;
- auto mesh_prop_velocity =
- const_cast<MeshLib::Mesh&>(mesh)
- .getProperties()
- .template createNewPropertyVector<double>(
- "velocity", MeshLib::MeshItemType::Cell, 3);
+ auto mesh_prop_velocity = MeshLib::getOrCreateMeshProperty<double>(
+ const_cast<MeshLib::Mesh&>(mesh),
+ "velocity", MeshLib::MeshItemType::Cell, 3);
mesh_prop_velocity->resize(mesh.getNumberOfElements() * 3);
_process_data.mesh_prop_velocity = mesh_prop_velocity;
if (!_vec_fracture_elements.empty())
{
- auto mesh_prop_levelset =
- const_cast<MeshLib::Mesh&>(mesh)
- .getProperties()
- .template createNewPropertyVector<double>(
- "levelset1", MeshLib::MeshItemType::Cell);
+ auto mesh_prop_levelset = MeshLib::getOrCreateMeshProperty<double>(
+ const_cast<MeshLib::Mesh&>(mesh),
+ "levelset1", MeshLib::MeshItemType::Cell, 1);
mesh_prop_levelset->resize(mesh.getNumberOfElements());
for (MeshLib::Element const* e : _mesh.getElements())
{
@@ -279,23 +266,20 @@ void HydroMechanicsProcess<GlobalDim>::initializeConcreteProcess(
(*mesh_prop_levelset)[e->getID()] = levelsets;
}
- auto mesh_prop_w_n = const_cast<MeshLib::Mesh&>(mesh)
- .getProperties()
- .template createNewPropertyVector<double>(
- "w_n", MeshLib::MeshItemType::Cell);
+ 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 = const_cast<MeshLib::Mesh&>(mesh)
- .getProperties()
- .template createNewPropertyVector<double>(
- "w_s", MeshLib::MeshItemType::Cell);
+ 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;
- auto mesh_prop_b = const_cast<MeshLib::Mesh&>(mesh)
- .getProperties()
- .template createNewPropertyVector<double>(
- "aperture", MeshLib::MeshItemType::Cell);
+ auto mesh_prop_b = MeshLib::getOrCreateMeshProperty<double>(
+ const_cast<MeshLib::Mesh&>(mesh),
+ "aperture", MeshLib::MeshItemType::Cell, 1);
mesh_prop_b->resize(mesh.getNumberOfElements());
auto mesh_prop_matid =
mesh.getProperties().getPropertyVector<int>("MaterialIDs");
@@ -312,39 +296,45 @@ void HydroMechanicsProcess<GlobalDim>::initializeConcreteProcess(
}
_process_data.mesh_prop_b = mesh_prop_b;
- auto mesh_prop_k_f = const_cast<MeshLib::Mesh&>(mesh)
- .getProperties()
- .template createNewPropertyVector<double>(
- "k_f", MeshLib::MeshItemType::Cell);
+
+ auto mesh_prop_k_f = MeshLib::getOrCreateMeshProperty<double>(
+ const_cast<MeshLib::Mesh&>(mesh),
+ "k_f", 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 =
- const_cast<MeshLib::Mesh&>(mesh)
- .getProperties()
- .template createNewPropertyVector<double>(
- "f_stress_s", MeshLib::MeshItemType::Cell);
+ 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 =
- const_cast<MeshLib::Mesh&>(mesh)
- .getProperties()
- .template createNewPropertyVector<double>(
- "f_stress_n", MeshLib::MeshItemType::Cell);
+ 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 =
- const_cast<MeshLib::Mesh&>(mesh)
- .getProperties()
- .template createNewPropertyVector<double>(
- "f_shear_failure", MeshLib::MeshItemType::Cell);
+ auto mesh_prop_fracture_shear_failure = MeshLib::getOrCreateMeshProperty<double>(
+ const_cast<MeshLib::Mesh&>(mesh),
+ "f_shear_failure", MeshLib::MeshItemType::Cell, 1);
mesh_prop_fracture_shear_failure->resize(mesh.getNumberOfElements());
_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;
}
}
@@ -357,6 +347,78 @@ void HydroMechanicsProcess<GlobalDim>::computeSecondaryVariableConcrete(
GlobalExecutor::executeMemberOnDereferenced(
&HydroMechanicsLocalAssemblerInterface::computeSecondaryVariable,
_local_assemblers, *_local_to_global_index_map, t, x, coupled_term);
+
+ // Copy displacement jumps in a solution vector to mesh property
+ // Remark: the copy is required because mesh properties for primary variables are set
+ // during output and are not ready yet when this function is called.
+ int g_variable_id = 0;
+ int g_global_component_offset = 0;
+ {
+ int global_component_offset_next = 0;
+ int global_component_offset = 0;
+ for (int variable_id = 0;
+ variable_id < static_cast<int>(this->getProcessVariables().size());
+ ++variable_id)
+ {
+ ProcessVariable& pv = this->getProcessVariables()[variable_id];
+ int const n_components = pv.getNumberOfComponents();
+ global_component_offset = global_component_offset_next;
+ global_component_offset_next += n_components;
+ if (pv.getName() != "displacement_jump1")
+ continue;
+
+ g_variable_id = variable_id;
+ g_global_component_offset = global_component_offset;
+ break;
+ }
+ }
+
+ MathLib::LinAlg::setLocalAccessibleVector(x);
+
+ ProcessVariable& pv_g = this->getProcessVariables()[g_variable_id];
+ auto& mesh_prop_g = pv_g.getOrCreateMeshProperty();
+ auto const num_comp = pv_g.getNumberOfComponents();
+ for (int component_id = 0; component_id < num_comp; ++component_id)
+ {
+ auto const& mesh_subsets =
+ _local_to_global_index_map->getMeshSubsets(g_variable_id,
+ component_id);
+ for (auto const& mesh_subset : mesh_subsets)
+ {
+ auto const mesh_id = mesh_subset->getMeshID();
+ for (auto const* node : mesh_subset->getNodes())
+ {
+ MeshLib::Location const l(
+ mesh_id, MeshLib::MeshItemType::Node, node->getID());
+
+ auto const global_component_id = g_global_component_offset + component_id;
+ mesh_prop_g[node->getID() * num_comp + component_id] =
+ x[global_component_id];
+ }
+ }
+ }
+
+ // compute nodal w and aperture
+ auto const& R = _process_data.fracture_property->R;
+ MeshLib::PropertyVector<double>& vec_w = *_process_data.mesh_prop_nodal_w;
+ MeshLib::PropertyVector<double>& vec_b = *_process_data.mesh_prop_nodal_b;
+
+ Eigen::VectorXd g(GlobalDim), w(GlobalDim);
+ for (MeshLib::Node const* node : _vec_fracture_nodes)
+ {
+ auto const node_id = node->getID();
+ g.setZero();
+ for (unsigned k=0; k<GlobalDim; k++)
+ g[k] = mesh_prop_g[node_id*GlobalDim + k];
+
+ w.noalias() = R * g;
+ for (unsigned k=0; k<GlobalDim; k++)
+ vec_w[node_id*GlobalDim + k] = w[k];
+
+ ProcessLib::SpatialPosition x;
+ x.setNodeID(node_id);
+ vec_b[node_id] = w[GlobalDim==2 ? 1 : 2] + (*_process_data.fracture_property->aperture0)(0,x)[0];
+ }
}
// ------------------------------------------------------------------------------------
diff --git a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcessData.h b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcessData.h
index f86382fe908cc4d5722ea5e3fbc250c59f6ad891..8d66a96936584e2d453913637503e0aa9e86d19b 100644
--- a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcessData.h
+++ b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcessData.h
@@ -140,6 +140,8 @@ struct HydroMechanicsProcessData
MeshLib::PropertyVector<double>* mesh_prop_fracture_stress_shear = 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;
};
} // namespace HydroMechanics