Skip to content
Snippets Groups Projects
Commit dbe8bbb3 authored by wenqing's avatar wenqing
Browse files

[HM] Added a member to HydroMechanicsLocalAssembler to assemble 

local matrices and vector for the staggered scheme
parent b02d6cfb
No related branches found
No related tags found
1 merge request!2016Staggered scheme for coupled processes with different orders of elements
Hide whitespace changes
Inline Side-by-side
ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h 0 → 100644
+ 248
− 0
View file @ dbe8bbb3
/**
* \copyright
* Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*
* \file HydroMechanicsFEM-impl.h
* Created on November 29, 2017, 2:03 PM
*/
#pragma once
#include "HydroMechanicsFEM.h"
#include "MaterialLib/SolidModels/KelvinVector.h"
#include "NumLib/Function/Interpolation.h"
namespace ProcessLib
{
namespace HydroMechanics
{
template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
typename IntegrationMethod, int DisplacementDim>
void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
ShapeFunctionPressure, IntegrationMethod,
DisplacementDim>::
assembleWithJacobianForPressureEquations(
const double t, const std::vector<double>& local_xdot,
const double /*dxdot_dx*/, const double /*dx_dx*/,
std::vector<double>& /*local_M_data*/,
std::vector<double>& /*local_K_data*/,
std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
const LocalCoupledSolutions& local_coupled_solutions)
{
auto const& local_p = local_coupled_solutions.local_coupled_xs[0];
auto const& local_u = local_coupled_solutions.local_coupled_xs[1];
assert(local_p.size() == pressure_size);
assert(local_u.size() == displacement_size);
auto const local_matrix_size = local_p.size();
auto local_rhs =
MathLib::createZeroedVector<typename ShapeMatricesTypeDisplacement::
template VectorType<pressure_size>>(
local_b_data, local_matrix_size);
SpatialPosition pos;
pos.setElementID(this->_element.getID());
auto p = Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
pressure_size> const>(local_p.data(), pressure_size);
auto u =
Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
displacement_size> const>(local_u.data(), displacement_size);
auto p_dot =
Eigen::Map<typename ShapeMatricesTypePressure::template VectorType<
pressure_size> const>(local_xdot.data(), pressure_size);
auto local_Jac = MathLib::createZeroedMatrix<
typename ShapeMatricesTypeDisplacement::template MatrixType<
pressure_size, pressure_size>>(local_Jac_data, pressure_size,
pressure_size);
typename ShapeMatricesTypePressure::NodalMatrixType mass =
ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
pressure_size);
typename ShapeMatricesTypePressure::NodalMatrixType laplace =
ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
pressure_size);
double const& dt = _process_data.dt;
SpatialPosition x_position;
x_position.setElementID(_element.getID());
int const n_integration_points = _integration_method.getNumberOfPoints();
for (int ip = 0; ip < n_integration_points; ip++)
{
x_position.setIntegrationPoint(ip);
auto const& w = _ip_data[ip].integration_weight;
auto const& N_u = _ip_data[ip].N_u;
auto const& dNdx_u = _ip_data[ip].dNdx_u;
auto const& N_p = _ip_data[ip].N_p;
auto const& dNdx_p = _ip_data[ip].dNdx_p;
double const S = _process_data.specific_storage(t, x_position)[0];
double const K_over_mu =
_process_data.intrinsic_permeability(t, x_position)[0] /
_process_data.fluid_viscosity(t, x_position)[0];
auto const alpha_b = _process_data.biot_coefficient(t, x_position)[0];
auto const rho_fr = _process_data.fluid_density(t, x_position)[0];
laplace.noalias() += dNdx_p.transpose() * K_over_mu * dNdx_p * w;
mass.noalias() += N_p.transpose() * S * N_p * w;
auto const& b = _process_data.specific_body_force;
local_rhs.noalias() += dNdx_p.transpose() * rho_fr * K_over_mu * b * w;
auto& eps = _ip_data[ip].eps;
auto const x_coord =
interpolateXCoordinate<ShapeFunctionDisplacement,
ShapeMatricesTypeDisplacement>(_element,
N_u);
auto const B =
LinearBMatrix::computeBMatrix<DisplacementDim,
ShapeFunctionDisplacement::NPOINTS,
typename BMatricesType::BMatrixType>(
dNdx_u, N_u, x_coord, _is_axially_symmetric);
eps.noalias() = B * u;
auto& eps_prev = _ip_data[ip].eps_prev;
const double dv_dt =
(Invariants::trace(eps) - Invariants::trace(eps_prev)) / dt;
local_rhs.noalias() -= alpha_b * dv_dt * N_p * w;
}
local_Jac.noalias() = laplace + mass / dt;
local_rhs.noalias() -= laplace * p + mass * p_dot;
}
template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
typename IntegrationMethod, int DisplacementDim>
void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
ShapeFunctionPressure, IntegrationMethod,
DisplacementDim>::
assembleWithJacobianForDeformationEquations(
const double t, const std::vector<double>& /*local_xdot*/,
const double /*dxdot_dx*/, const double /*dx_dx*/,
std::vector<double>& /*local_M_data*/,
std::vector<double>& /*local_K_data*/,
std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
const LocalCoupledSolutions& local_coupled_solutions)
{
auto const& local_p = local_coupled_solutions.local_coupled_xs[0];
auto const& local_u = local_coupled_solutions.local_coupled_xs[1];
assert(local_p.size() == pressure_size);
assert(local_u.size() == displacement_size);
auto u =
Eigen::Map<typename ShapeMatricesTypeDisplacement::template VectorType<
displacement_size> const>(local_u.data(), displacement_size);
auto local_Jac = MathLib::createZeroedMatrix<
typename ShapeMatricesTypeDisplacement::template MatrixType<
displacement_size, displacement_size>>(
local_Jac_data, displacement_size, displacement_size);
auto local_rhs =
MathLib::createZeroedVector<typename ShapeMatricesTypeDisplacement::
template VectorType<displacement_size>>(
local_b_data, displacement_size);
double const& dt = _process_data.dt;
SpatialPosition x_position;
x_position.setElementID(_element.getID());
int const n_integration_points = _integration_method.getNumberOfPoints();
for (int ip = 0; ip < n_integration_points; ip++)
{
x_position.setIntegrationPoint(ip);
auto const& w = _ip_data[ip].integration_weight;
auto const& N_u_op = _ip_data[ip].N_u_op;
auto const& N_u = _ip_data[ip].N_u;
auto const& dNdx_u = _ip_data[ip].dNdx_u;
auto const& N_p = _ip_data[ip].N_p;
auto const x_coord =
interpolateXCoordinate<ShapeFunctionDisplacement,
ShapeMatricesTypeDisplacement>(_element,
N_u);
auto const B =
LinearBMatrix::computeBMatrix<DisplacementDim,
ShapeFunctionDisplacement::NPOINTS,
typename BMatricesType::BMatrixType>(
dNdx_u, N_u, x_coord, _is_axially_symmetric);
auto& eps = _ip_data[ip].eps;
auto const& sigma_eff = _ip_data[ip].sigma_eff;
auto const alpha = _process_data.biot_coefficient(t, x_position)[0];
auto const rho_sr = _process_data.solid_density(t, x_position)[0];
auto const rho_fr = _process_data.fluid_density(t, x_position)[0];
auto const porosity = _process_data.porosity(t, x_position)[0];
auto const& b = _process_data.specific_body_force;
auto const& identity2 = MaterialLib::SolidModels::Invariants<
KelvinVectorDimensions<DisplacementDim>::value>::identity2;
eps.noalias() = B * u;
auto C = _ip_data[ip].updateConstitutiveRelation(t, x_position, dt, u);
local_Jac
.template block<displacement_size, displacement_size>(
displacement_index, displacement_index)
.noalias() += B.transpose() * C * B * w;
double p_at_xi = 0.;
NumLib::shapeFunctionInterpolate(local_p, N_p, p_at_xi);
double const rho = rho_sr * (1 - porosity) + porosity * rho_fr;
local_rhs.template segment<displacement_size>(displacement_index)
.noalias() -=
(B.transpose() * (sigma_eff - alpha * identity2 * p_at_xi) -
N_u_op.transpose() * rho * b) *
w;
}
}
template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
typename IntegrationMethod, int DisplacementDim>
void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
ShapeFunctionPressure, IntegrationMethod,
DisplacementDim>::
assembleWithJacobianAndCoupling(
const double t,
const std::vector<double>& local_xdot,
const double dxdot_dx,
const double dx_dx,
std::vector<double>& local_M_data,
std::vector<double>& local_K_data,
std::vector<double>& local_b_data,
std::vector<double>& local_Jac_data,
const LocalCoupledSolutions& local_coupled_solutions)
{
// For the equations with pressure
if (local_coupled_solutions.process_id == 0)
{
assembleWithJacobianForPressureEquations(
t, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
local_b_data, local_Jac_data, local_coupled_solutions);
return;
}
// For the equations with deformation
assembleWithJacobianForDeformationEquations(
t, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
local_b_data, local_Jac_data, local_coupled_solutions);
}
} // namespace HydroMechanics
} // namespace ProcessLib
ProcessLib/HydroMechanics/HydroMechanicsFEM.h
+ 29
− 1
View file @ dbe8bbb3
...@@ -110,6 +110,10 @@ public: ...@@ -110,6 +110,10 @@ public:
using ShapeMatricesTypePressure = using ShapeMatricesTypePressure =
ShapeMatrixPolicyType<ShapeFunctionPressure, DisplacementDim>; ShapeMatrixPolicyType<ShapeFunctionPressure, DisplacementDim>;
static int const KelvinVectorSize =
ProcessLib::KelvinVectorDimensions<DisplacementDim>::value;
using Invariants = MaterialLib::SolidModels::Invariants<KelvinVectorSize>;
HydroMechanicsLocalAssembler(HydroMechanicsLocalAssembler const&) = delete; HydroMechanicsLocalAssembler(HydroMechanicsLocalAssembler const&) = delete;
HydroMechanicsLocalAssembler(HydroMechanicsLocalAssembler&&) = delete; HydroMechanicsLocalAssembler(HydroMechanicsLocalAssembler&&) = delete;
...@@ -352,6 +356,13 @@ public: ...@@ -352,6 +356,13 @@ public:
.noalias() += Kup * p; .noalias() += Kup * p;
} }
void assembleWithJacobianAndCoupling(
double const t, std::vector<double> const& local_xdot,
const double dxdot_dx, const double dx_dx,
std::vector<double>& local_M_data, std::vector<double>& local_K_data,
std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
LocalCoupledSolutions const& local_coupled_solutions) override;
void preTimestepConcrete(std::vector<double> const& /*local_x*/, void preTimestepConcrete(std::vector<double> const& /*local_x*/,
double const /*t*/, double const /*t*/,
double const /*delta_t*/) override double const /*delta_t*/) override
...@@ -516,7 +527,8 @@ public: ...@@ -516,7 +527,8 @@ public:
SpatialPosition x_position; SpatialPosition x_position;
x_position.setElementID(_element.getID()); x_position.setElementID(_element.getID());
for (unsigned ip = 0; ip < n_integration_points; ip++) { for (unsigned ip = 0; ip < n_integration_points; ip++)
{
x_position.setIntegrationPoint(ip); x_position.setIntegrationPoint(ip);
double const K_over_mu = double const K_over_mu =
_process_data.intrinsic_permeability(t, x_position)[0] / _process_data.intrinsic_permeability(t, x_position)[0] /
...@@ -566,6 +578,20 @@ private: ...@@ -566,6 +578,20 @@ private:
return cache; return cache;
} }
void assembleWithJacobianForDeformationEquations(
double const t, std::vector<double> const& local_xdot,
const double dxdot_dx, const double dx_dx,
std::vector<double>& local_M_data, std::vector<double>& local_K_data,
std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
LocalCoupledSolutions const& local_coupled_solutions);
void assembleWithJacobianForPressureEquations(
double const t, std::vector<double> const& local_xdot,
const double dxdot_dx, const double dx_dx,
std::vector<double>& local_M_data, std::vector<double>& local_K_data,
std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
LocalCoupledSolutions const& local_coupled_solutions);
private: private:
HydroMechanicsProcessData<DisplacementDim>& _process_data; HydroMechanicsProcessData<DisplacementDim>& _process_data;
...@@ -595,3 +621,5 @@ private: ...@@ -595,3 +621,5 @@ private:
} // namespace HydroMechanics } // namespace HydroMechanics
} // namespace ProcessLib } // namespace ProcessLib
#include "HydroMechanicsFEM-impl.h"
\ No newline at end of file
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