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

Merge pull request #1859 from norihiro-w/lie-add-nodal-w 

[PCS/LIE/HM] add output of nodal w and aperture
parents 54874f44 3b057f1c
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp
+ 130
− 68
View file @ b7eaf3d7
...@@ -207,66 +207,53 @@ void HydroMechanicsProcess<GlobalDim>::initializeConcreteProcess( ...@@ -207,66 +207,53 @@ void HydroMechanicsProcess<GlobalDim>::initializeConcreteProcess(
_local_assemblers, mesh.isAxiallySymmetric(), integration_order, _local_assemblers, mesh.isAxiallySymmetric(), integration_order,
_process_data); _process_data);
auto mesh_prop_sigma_xx = const_cast<MeshLib::Mesh&>(mesh) auto mesh_prop_sigma_xx = MeshLib::getOrCreateMeshProperty<double>(
.getProperties() const_cast<MeshLib::Mesh&>(mesh),
.template createNewPropertyVector<double>( "stress_xx", MeshLib::MeshItemType::Cell, 1);
"stress_xx", MeshLib::MeshItemType::Cell);
mesh_prop_sigma_xx->resize(mesh.getNumberOfElements()); mesh_prop_sigma_xx->resize(mesh.getNumberOfElements());
_process_data.mesh_prop_stress_xx = mesh_prop_sigma_xx; _process_data.mesh_prop_stress_xx = mesh_prop_sigma_xx;
auto mesh_prop_sigma_yy = const_cast<MeshLib::Mesh&>(mesh) auto mesh_prop_sigma_yy = MeshLib::getOrCreateMeshProperty<double>(
.getProperties() const_cast<MeshLib::Mesh&>(mesh),
.template createNewPropertyVector<double>( "stress_yy", MeshLib::MeshItemType::Cell, 1);
"stress_yy", MeshLib::MeshItemType::Cell);
mesh_prop_sigma_yy->resize(mesh.getNumberOfElements()); mesh_prop_sigma_yy->resize(mesh.getNumberOfElements());
_process_data.mesh_prop_stress_yy = mesh_prop_sigma_yy; _process_data.mesh_prop_stress_yy = mesh_prop_sigma_yy;
auto mesh_prop_sigma_xy = const_cast<MeshLib::Mesh&>(mesh) auto mesh_prop_sigma_xy = MeshLib::getOrCreateMeshProperty<double>(
.getProperties() const_cast<MeshLib::Mesh&>(mesh),
.template createNewPropertyVector<double>( "stress_xy", MeshLib::MeshItemType::Cell, 1);
"stress_xy", MeshLib::MeshItemType::Cell);
mesh_prop_sigma_xy->resize(mesh.getNumberOfElements()); mesh_prop_sigma_xy->resize(mesh.getNumberOfElements());
_process_data.mesh_prop_stress_xy = mesh_prop_sigma_xy; _process_data.mesh_prop_stress_xy = mesh_prop_sigma_xy;
auto mesh_prop_epsilon_xx = auto mesh_prop_epsilon_xx = MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh) const_cast<MeshLib::Mesh&>(mesh),
.getProperties() "strain_xx", MeshLib::MeshItemType::Cell, 1);
.template createNewPropertyVector<double>(
"strain_xx", MeshLib::MeshItemType::Cell);
mesh_prop_epsilon_xx->resize(mesh.getNumberOfElements()); mesh_prop_epsilon_xx->resize(mesh.getNumberOfElements());
_process_data.mesh_prop_strain_xx = mesh_prop_epsilon_xx; _process_data.mesh_prop_strain_xx = mesh_prop_epsilon_xx;
auto mesh_prop_epsilon_yy = auto mesh_prop_epsilon_yy = MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh) const_cast<MeshLib::Mesh&>(mesh),
.getProperties() "strain_yy", MeshLib::MeshItemType::Cell, 1);
.template createNewPropertyVector<double>(
"strain_yy", MeshLib::MeshItemType::Cell);
mesh_prop_epsilon_yy->resize(mesh.getNumberOfElements()); mesh_prop_epsilon_yy->resize(mesh.getNumberOfElements());
_process_data.mesh_prop_strain_yy = mesh_prop_epsilon_yy; _process_data.mesh_prop_strain_yy = mesh_prop_epsilon_yy;
auto mesh_prop_epsilon_xy = auto mesh_prop_epsilon_xy = MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh) const_cast<MeshLib::Mesh&>(mesh),
.getProperties() "strain_xy", MeshLib::MeshItemType::Cell, 1);
.template createNewPropertyVector<double>(
"strain_xy", MeshLib::MeshItemType::Cell);
mesh_prop_epsilon_xy->resize(mesh.getNumberOfElements()); mesh_prop_epsilon_xy->resize(mesh.getNumberOfElements());
_process_data.mesh_prop_strain_xy = mesh_prop_epsilon_xy; _process_data.mesh_prop_strain_xy = mesh_prop_epsilon_xy;
auto mesh_prop_velocity = auto mesh_prop_velocity = MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh) const_cast<MeshLib::Mesh&>(mesh),
.getProperties() "velocity", MeshLib::MeshItemType::Cell, 3);
.template createNewPropertyVector<double>(
"velocity", MeshLib::MeshItemType::Cell, 3);
mesh_prop_velocity->resize(mesh.getNumberOfElements() * 3); mesh_prop_velocity->resize(mesh.getNumberOfElements() * 3);
_process_data.mesh_prop_velocity = mesh_prop_velocity; _process_data.mesh_prop_velocity = mesh_prop_velocity;
if (!_vec_fracture_elements.empty()) if (!_vec_fracture_elements.empty())
{ {
auto mesh_prop_levelset = auto mesh_prop_levelset = MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh) const_cast<MeshLib::Mesh&>(mesh),
.getProperties() "levelset1", MeshLib::MeshItemType::Cell, 1);
.template createNewPropertyVector<double>(
"levelset1", MeshLib::MeshItemType::Cell);
mesh_prop_levelset->resize(mesh.getNumberOfElements()); mesh_prop_levelset->resize(mesh.getNumberOfElements());
for (MeshLib::Element const* e : _mesh.getElements()) for (MeshLib::Element const* e : _mesh.getElements())
{ {
...@@ -279,23 +266,20 @@ void HydroMechanicsProcess<GlobalDim>::initializeConcreteProcess( ...@@ -279,23 +266,20 @@ void HydroMechanicsProcess<GlobalDim>::initializeConcreteProcess(
(*mesh_prop_levelset)[e->getID()] = levelsets; (*mesh_prop_levelset)[e->getID()] = levelsets;
} }
auto mesh_prop_w_n = const_cast<MeshLib::Mesh&>(mesh) auto mesh_prop_w_n = MeshLib::getOrCreateMeshProperty<double>(
.getProperties() const_cast<MeshLib::Mesh&>(mesh),
.template createNewPropertyVector<double>( "w_n", MeshLib::MeshItemType::Cell, 1);
"w_n", MeshLib::MeshItemType::Cell);
mesh_prop_w_n->resize(mesh.getNumberOfElements()); mesh_prop_w_n->resize(mesh.getNumberOfElements());
auto mesh_prop_w_s = const_cast<MeshLib::Mesh&>(mesh) auto mesh_prop_w_s = MeshLib::getOrCreateMeshProperty<double>(
.getProperties() const_cast<MeshLib::Mesh&>(mesh),
.template createNewPropertyVector<double>( "w_s", MeshLib::MeshItemType::Cell, 1);
"w_s", MeshLib::MeshItemType::Cell);
mesh_prop_w_s->resize(mesh.getNumberOfElements()); mesh_prop_w_s->resize(mesh.getNumberOfElements());
_process_data.mesh_prop_w_n = mesh_prop_w_n; _process_data.mesh_prop_w_n = mesh_prop_w_n;
_process_data.mesh_prop_w_s = mesh_prop_w_s; _process_data.mesh_prop_w_s = mesh_prop_w_s;
auto mesh_prop_b = const_cast<MeshLib::Mesh&>(mesh) auto mesh_prop_b = MeshLib::getOrCreateMeshProperty<double>(
.getProperties() const_cast<MeshLib::Mesh&>(mesh),
.template createNewPropertyVector<double>( "aperture", MeshLib::MeshItemType::Cell, 1);
"aperture", MeshLib::MeshItemType::Cell);
mesh_prop_b->resize(mesh.getNumberOfElements()); mesh_prop_b->resize(mesh.getNumberOfElements());
auto mesh_prop_matid = auto mesh_prop_matid =
mesh.getProperties().getPropertyVector<int>("MaterialIDs"); mesh.getProperties().getPropertyVector<int>("MaterialIDs");
...@@ -312,39 +296,45 @@ void HydroMechanicsProcess<GlobalDim>::initializeConcreteProcess( ...@@ -312,39 +296,45 @@ void HydroMechanicsProcess<GlobalDim>::initializeConcreteProcess(
} }
_process_data.mesh_prop_b = mesh_prop_b; _process_data.mesh_prop_b = mesh_prop_b;
auto mesh_prop_k_f = const_cast<MeshLib::Mesh&>(mesh)
.getProperties() auto mesh_prop_k_f = MeshLib::getOrCreateMeshProperty<double>(
.template createNewPropertyVector<double>( const_cast<MeshLib::Mesh&>(mesh),
"k_f", MeshLib::MeshItemType::Cell); "k_f", MeshLib::MeshItemType::Cell, 1);
mesh_prop_k_f->resize(mesh.getNumberOfElements()); mesh_prop_k_f->resize(mesh.getNumberOfElements());
_process_data.mesh_prop_k_f = mesh_prop_k_f; _process_data.mesh_prop_k_f = mesh_prop_k_f;
auto mesh_prop_fracture_stress_shear = auto mesh_prop_fracture_stress_shear = MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh) const_cast<MeshLib::Mesh&>(mesh),
.getProperties() "f_stress_s", MeshLib::MeshItemType::Cell, 1);
.template createNewPropertyVector<double>(
"f_stress_s", MeshLib::MeshItemType::Cell);
mesh_prop_fracture_stress_shear->resize(mesh.getNumberOfElements()); mesh_prop_fracture_stress_shear->resize(mesh.getNumberOfElements());
_process_data.mesh_prop_fracture_stress_shear = _process_data.mesh_prop_fracture_stress_shear =
mesh_prop_fracture_stress_shear; mesh_prop_fracture_stress_shear;
auto mesh_prop_fracture_stress_normal = auto mesh_prop_fracture_stress_normal = MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh) const_cast<MeshLib::Mesh&>(mesh),
.getProperties() "f_stress_n", MeshLib::MeshItemType::Cell, 1);
.template createNewPropertyVector<double>(
"f_stress_n", MeshLib::MeshItemType::Cell);
mesh_prop_fracture_stress_normal->resize(mesh.getNumberOfElements()); mesh_prop_fracture_stress_normal->resize(mesh.getNumberOfElements());
_process_data.mesh_prop_fracture_stress_normal = _process_data.mesh_prop_fracture_stress_normal =
mesh_prop_fracture_stress_normal; mesh_prop_fracture_stress_normal;
auto mesh_prop_fracture_shear_failure = auto mesh_prop_fracture_shear_failure = MeshLib::getOrCreateMeshProperty<double>(
const_cast<MeshLib::Mesh&>(mesh) const_cast<MeshLib::Mesh&>(mesh),
.getProperties() "f_shear_failure", MeshLib::MeshItemType::Cell, 1);
.template createNewPropertyVector<double>(
"f_shear_failure", MeshLib::MeshItemType::Cell);
mesh_prop_fracture_shear_failure->resize(mesh.getNumberOfElements()); mesh_prop_fracture_shear_failure->resize(mesh.getNumberOfElements());
_process_data.mesh_prop_fracture_shear_failure = _process_data.mesh_prop_fracture_shear_failure =
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( ...@@ -357,6 +347,78 @@ void HydroMechanicsProcess<GlobalDim>::computeSecondaryVariableConcrete(
GlobalExecutor::executeMemberOnDereferenced( GlobalExecutor::executeMemberOnDereferenced(
&HydroMechanicsLocalAssemblerInterface::computeSecondaryVariable, &HydroMechanicsLocalAssemblerInterface::computeSecondaryVariable,
_local_assemblers, *_local_to_global_index_map, t, x, coupled_term); _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];
}
} }
// ------------------------------------------------------------------------------------ // ------------------------------------------------------------------------------------
......
ProcessLib/LIE/HydroMechanics/HydroMechanicsProcessData.h
+ 2
− 0
View file @ b7eaf3d7
...@@ -140,6 +140,8 @@ struct HydroMechanicsProcessData ...@@ -140,6 +140,8 @@ struct HydroMechanicsProcessData
MeshLib::PropertyVector<double>* mesh_prop_fracture_stress_shear = nullptr; MeshLib::PropertyVector<double>* mesh_prop_fracture_stress_shear = nullptr;
MeshLib::PropertyVector<double>* mesh_prop_fracture_stress_normal = nullptr; MeshLib::PropertyVector<double>* mesh_prop_fracture_stress_normal = nullptr;
MeshLib::PropertyVector<double>* mesh_prop_fracture_shear_failure = 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 } // namespace HydroMechanics
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment