diff --git a/ProcessLib/PhaseField/CreatePhaseFieldProcess.cpp b/ProcessLib/PhaseField/CreatePhaseFieldProcess.cpp
index 8116905574936e1471daebded5a895bc4545687b..515427e6874eb6b69c7814fb2e1f1e96b248dc65 100644
--- a/ProcessLib/PhaseField/CreatePhaseFieldProcess.cpp
+++ b/ProcessLib/PhaseField/CreatePhaseFieldProcess.cpp
@@ -22,7 +22,6 @@ namespace ProcessLib
{
namespace PhaseField
{
-
template <int DisplacementDim>
std::unique_ptr<Process> createPhaseFieldProcess(
MeshLib::Mesh& mesh,
@@ -56,9 +55,9 @@ std::unique_ptr<Process> createPhaseFieldProcess(
auto per_process_variables = findProcessVariables(
variables, pv_config,
{//! \ogs_file_param_special{prj__processes__process__PHASE_FIELD__process_variables__phasefield}
- "phasefield",
+ "phasefield",
//! \ogs_file_param_special{prj__processes__process__PHASE_FIELD__process_variables__displacement}
- "displacement"});
+ "displacement"});
variable_ph = &per_process_variables[0].get();
variable_u = &per_process_variables[1].get();
process_variables.push_back(std::move(per_process_variables));
@@ -66,14 +65,14 @@ std::unique_ptr<Process> createPhaseFieldProcess(
else // staggered scheme.
{
using namespace std::string_literals;
- for (auto const& variable_name : {"phasefield"s, "displacement"s})
+ for (auto const& variable_name : {"displacement"s, "phasefield"s})
{
auto per_process_variables =
findProcessVariables(variables, pv_config, {variable_name});
process_variables.push_back(std::move(per_process_variables));
}
- variable_ph = &process_variables[0][0].get();
- variable_u = &process_variables[1][0].get();
+ variable_u = &process_variables[0][0].get();
+ variable_ph = &process_variables[1][0].get();
}
DBUG("Associate displacement with process variable \'%s\'.",
@@ -94,7 +93,7 @@ std::unique_ptr<Process> createPhaseFieldProcess(
if (variable_ph->getNumberOfComponents() != 1)
{
OGS_FATAL(
- "Pressure process variable '%s' is not a scalar variable but has "
+ "Phasefield process variable '%s' is not a scalar variable but has "
"%d components.",
variable_ph->getName().c_str(),
variable_ph->getNumberOfComponents());
@@ -118,11 +117,13 @@ std::unique_ptr<Process> createPhaseFieldProcess(
material = nullptr;
if (type == "LinearElasticIsotropic")
{
- auto elastic_model = MaterialLib::Solids::createLinearElasticIsotropic<
- DisplacementDim>(parameters, constitutive_relation_config);
- material =
- std::make_unique<MaterialLib::Solids::LinearElasticIsotropicPhaseField<
- DisplacementDim>>(std::move(elastic_model->getMaterialProperties()));
+ auto elastic_model =
+ MaterialLib::Solids::createLinearElasticIsotropic<DisplacementDim>(
+ parameters, constitutive_relation_config);
+ material = std::make_unique<
+ MaterialLib::Solids::LinearElasticIsotropicPhaseField<
+ DisplacementDim>>(
+ std::move(elastic_model->getMaterialProperties()));
}
else
{
@@ -205,10 +206,10 @@ std::unique_ptr<Process> createPhaseFieldProcess(
named_function_caller);
return std::make_unique<PhaseFieldProcess<DisplacementDim>>(
- mesh, std::move(jacobian_assembler), parameters, integration_order,
- std::move(process_variables), std::move(process_data),
- std::move(secondary_variables), std::move(named_function_caller),
- use_monolithic_scheme);
+ mesh, std::move(jacobian_assembler), parameters, integration_order,
+ std::move(process_variables), std::move(process_data),
+ std::move(secondary_variables), std::move(named_function_caller),
+ use_monolithic_scheme);
}
template std::unique_ptr<Process> createPhaseFieldProcess<2>(
diff --git a/ProcessLib/PhaseField/PhaseFieldFEM-impl.h b/ProcessLib/PhaseField/PhaseFieldFEM-impl.h
index 812a4b36b1e6c5d554529354cec95ef3e29a68aa..5d3433ccc1a9a803f5822986f899f4089755ffac 100644
--- a/ProcessLib/PhaseField/PhaseFieldFEM-impl.h
+++ b/ProcessLib/PhaseField/PhaseFieldFEM-impl.h
@@ -207,7 +207,7 @@ void PhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
LocalCoupledSolutions const& local_coupled_solutions)
{
// For the equations with phase field.
- if (local_coupled_solutions.process_id == 0)
+ if (local_coupled_solutions.process_id == 1)
{
assembleWithJacobianPhaseFiledEquations(
t, local_xdot, dxdot_dx, dx_dx, local_M_data, local_K_data,
@@ -230,12 +230,12 @@ void PhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
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) const
+ LocalCoupledSolutions const& local_coupled_solutions)
{
- auto const& local_d = local_coupled_solutions.local_coupled_xs[0];
- auto const& local_u = local_coupled_solutions.local_coupled_xs[1];
- assert(local_d.size() == phasefield_size);
+ auto const& local_u = local_coupled_solutions.local_coupled_xs[0];
+ auto const& local_d = local_coupled_solutions.local_coupled_xs[1];
assert(local_u.size() == displacement_size);
+ assert(local_d.size() == phasefield_size);
auto const local_matrix_size = local_u.size();
auto d = Eigen::Map<
@@ -245,26 +245,15 @@ void PhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
auto u = Eigen::Map<typename ShapeMatricesType::template VectorType<
displacement_size> const>(local_u.data(), displacement_size);
- // May not needed: auto d_dot = Eigen::Map<
- // typename ShapeMatricesType::template VectorType<phasefield_size>
- // const>(
- // local_xdot.data(), phasefield_size);
-
auto local_Jac = MathLib::createZeroedMatrix<JacobianMatrix>(
local_Jac_data, local_matrix_size, local_matrix_size);
auto local_rhs =
MathLib::createZeroedVector<RhsVector>(local_b_data, local_matrix_size);
- typename ShapeMatricesType::template MatrixType<displacement_size,
- phasefield_size>
- Kud;
- Kud.setZero(displacement_size, phasefield_size);
-
- double const& dt = _process_data.dt;
-
SpatialPosition x_position;
x_position.setElementID(_element.getID());
+ double const& dt = _process_data.dt;
int const n_integration_points = _integration_method.getNumberOfPoints();
for (int ip = 0; ip < n_integration_points; ip++)
@@ -273,6 +262,52 @@ void PhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
auto const& w = _ip_data[ip].integration_weight;
auto const& N = _ip_data[ip].N;
auto const& dNdx = _ip_data[ip].dNdx;
+
+ auto const x_coord =
+ interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(_element,
+ N);
+
+ auto const& B =
+ LinearBMatrix::computeBMatrix<DisplacementDim,
+ ShapeFunction::NPOINTS,
+ typename BMatricesType::BMatrixType>(
+ dNdx, N, x_coord, _is_axially_symmetric);
+
+ auto& eps = _ip_data[ip].eps;
+ eps.noalias() = B * u;
+ double const k = _process_data.residual_stiffness(t, x_position)[0];
+ double const d_ip = N.dot(d);
+ double const degradation = d_ip * d_ip * (1 - k) + k;
+ _ip_data[ip].updateConstitutiveRelation(t, x_position, dt, u,
+ degradation);
+
+ auto& sigma_real = _ip_data[ip].sigma_real;
+ auto const& C_tensile = _ip_data[ip].C_tensile;
+ auto const& C_compressive = _ip_data[ip].C_compressive;
+
+ typename ShapeMatricesType::template MatrixType<DisplacementDim,
+ displacement_size>
+ N_u = ShapeMatricesType::template MatrixType<
+ DisplacementDim, displacement_size>::Zero(DisplacementDim,
+ displacement_size);
+
+ for (int i = 0; i < DisplacementDim; ++i)
+ N_u.template block<1, displacement_size / DisplacementDim>(
+ i, i * displacement_size / DisplacementDim)
+ .noalias() = N;
+
+ auto const rho_sr = _process_data.solid_density(t, x_position)[0];
+ auto const& b = _process_data.specific_body_force;
+
+ auto pressure = _ip_data[ip].pressure;
+
+ local_rhs.noalias() -=
+ (B.transpose() * sigma_real - N_u.transpose() * rho_sr * b -
+ pressure * N_u.transpose() * dNdx * d) *
+ w;
+
+ local_Jac.noalias() +=
+ B.transpose() * (degradation * C_tensile + C_compressive) * B * w;
}
}
@@ -285,12 +320,12 @@ void PhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
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) const
+ LocalCoupledSolutions const& local_coupled_solutions)
{
- auto const& local_d = local_coupled_solutions.local_coupled_xs[0];
- auto const& local_u = local_coupled_solutions.local_coupled_xs[1];
- assert(local_d.size() == phasefield_size);
+ auto const& local_u = local_coupled_solutions.local_coupled_xs[0];
+ auto const& local_d = local_coupled_solutions.local_coupled_xs[1];
assert(local_u.size() == displacement_size);
+ assert(local_d.size() == phasefield_size);
auto const local_matrix_size = local_d.size();
auto d = Eigen::Map<
@@ -300,29 +335,12 @@ void PhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
auto u = Eigen::Map<typename ShapeMatricesType::template VectorType<
displacement_size> const>(local_u.data(), displacement_size);
- auto d_dot = Eigen::Map<
- typename ShapeMatricesType::template VectorType<phasefield_size> const>(
- local_xdot.data(), phasefield_size);
-
auto local_Jac = MathLib::createZeroedMatrix<JacobianMatrix>(
local_Jac_data, local_matrix_size, local_matrix_size);
auto local_rhs =
MathLib::createZeroedVector<RhsVector>(local_b_data, local_matrix_size);
- typename ShapeMatricesType::template MatrixType<phasefield_size,
- displacement_size>
- Kdu;
- Kdu.setZero(phasefield_size, displacement_size);
-
- typename ShapeMatricesType::NodalMatrixType Kdd;
- Kdd.setZero(phasefield_size, phasefield_size);
-
- typename ShapeMatricesType::NodalMatrixType Ddd;
- Ddd.setZero(phasefield_size, phasefield_size);
-
- double const& dt = _process_data.dt;
-
SpatialPosition x_position;
x_position.setElementID(_element.getID());
@@ -333,8 +351,40 @@ void PhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
auto const& w = _ip_data[ip].integration_weight;
auto const& N = _ip_data[ip].N;
auto const& dNdx = _ip_data[ip].dNdx;
+
+ double const gc = _process_data.crack_resistance(t, x_position)[0];
+ double const ls = _process_data.crack_length_scale(t, x_position)[0];
+ auto const& strain_energy_tensile = _ip_data[ip].strain_energy_tensile;
+
+ auto& ip_data = _ip_data[ip];
+ ip_data.strain_energy_tensile = strain_energy_tensile;
+
+ typename ShapeMatricesType::template MatrixType<DisplacementDim,
+ displacement_size>
+ N_u = ShapeMatricesType::template MatrixType<
+ DisplacementDim, displacement_size>::Zero(DisplacementDim,
+ displacement_size);
+
+ for (int i = 0; i < DisplacementDim; ++i)
+ N_u.template block<1, displacement_size / DisplacementDim>(
+ i, i * displacement_size / DisplacementDim)
+ .noalias() = N;
+
+ auto pressure = _ip_data[ip].pressure;
+
+ local_Jac.noalias() +=
+ (2 * N.transpose() * N * strain_energy_tensile +
+ gc * (N.transpose() * N / ls + dNdx.transpose() * dNdx * ls)) *
+ w;
+
+ local_rhs.noalias() -=
+ (N.transpose() * N * d * 2 * strain_energy_tensile +
+ gc * ((N.transpose() * N / ls + dNdx.transpose() * dNdx * ls) * d -
+ N.transpose() / ls) -
+ pressure * dNdx.transpose() * N_u * u) *
+ w;
}
}
} // namespace PhaseField
-} // namespace ProcessLib
\ No newline at end of file
+} // namespace ProcessLib
diff --git a/ProcessLib/PhaseField/PhaseFieldFEM.h b/ProcessLib/PhaseField/PhaseFieldFEM.h
index 7486f72abf1b544f0d71b74995ee56cd3712501d..ea47df87a675d864af14d9b7b018fb35747f247f 100644
--- a/ProcessLib/PhaseField/PhaseFieldFEM.h
+++ b/ProcessLib/PhaseField/PhaseFieldFEM.h
@@ -57,6 +57,7 @@ struct IntegrationPointData final
double integration_weight;
double history_variable;
double history_variable_prev;
+ double pressure;
void pushBackState()
{
@@ -165,6 +166,11 @@ public:
ip_data.history_variable_prev =
_process_data.history_field(0, x_position)[0];
ip_data.sigma_real.setZero(kelvin_vector_size);
+ ip_data.strain_energy_tensile = 0.0;
+
+ /// pressure is initialized with unity (=1.0) as the problem is
+ /// linearly scalable in post-process.
+ ip_data.pressure = 1.0;
ip_data.N = shape_matrices[ip].N;
ip_data.dNdx = shape_matrices[ip].dNdx;
@@ -373,14 +379,14 @@ private:
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) const;
+ LocalCoupledSolutions const& local_coupled_solutions);
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) const;
+ LocalCoupledSolutions const& local_coupled_solutions);
PhaseFieldProcessData<DisplacementDim>& _process_data;
diff --git a/ProcessLib/PhaseField/PhaseFieldProcess-impl.h b/ProcessLib/PhaseField/PhaseFieldProcess-impl.h
index eaffca16539e06a7e5d27af1b5942b09c6412bee..a6d92bfab2c1c1796c93bf34603f3e2e6369e94d 100644
--- a/ProcessLib/PhaseField/PhaseFieldProcess-impl.h
+++ b/ProcessLib/PhaseField/PhaseFieldProcess-impl.h
@@ -55,7 +55,7 @@ PhaseFieldProcess<DisplacementDim>::getMatrixSpecifications(
{
// For the monolithic scheme or the M process (deformation) in the staggered
// scheme.
- if (_use_monolithic_scheme || process_id == 1)
+ if (_use_monolithic_scheme || process_id == 0)
{
auto const& l = *_local_to_global_index_map;
return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
@@ -74,12 +74,12 @@ PhaseFieldProcess<DisplacementDim>::getDOFTable(const int process_id) const
{
// If monolithic scheme is used or the equation of deformation is solved in
// the staggered scheme.
- if (_use_monolithic_scheme || process_id == 1)
+ if (_use_monolithic_scheme || process_id == 0)
{
return *_local_to_global_index_map;
}
- // For the equation of pressure
+ // For the equation of phasefield
return *_local_to_global_index_map_single_component;
}
@@ -125,7 +125,7 @@ void PhaseFieldProcess<DisplacementDim>::constructDofTable()
else
{
// For displacement equation.
- const int process_id = 1;
+ const int process_id = 0;
std::vector<MeshLib::MeshSubsets> all_mesh_subsets;
std::generate_n(
std::back_inserter(all_mesh_subsets),
@@ -230,22 +230,21 @@ void PhaseFieldProcess<DisplacementDim>::initializeBoundaryConditions()
{
if (_use_monolithic_scheme)
{
- const int process_id_of_pf = 0;
+ const int process_id_of_phasefieldmechanics = 0;
initializeProcessBoundaryConditionsAndSourceTerms(
- *_local_to_global_index_map, process_id_of_pf);
+ *_local_to_global_index_map, process_id_of_phasefieldmechanics);
return;
}
// Staggered scheme:
- // for the phase field
- const int process_id_of_pf = 0;
- initializeProcessBoundaryConditionsAndSourceTerms(
- *_local_to_global_index_map_single_component, process_id_of_pf);
-
// for the equations of deformation.
- const int process_id_of_u = 1;
+ const int mechanical_process_id = 0;
initializeProcessBoundaryConditionsAndSourceTerms(
- *_local_to_global_index_map, process_id_of_u);
+ *_local_to_global_index_map, mechanical_process_id);
+ // for the phase field
+ const int phasefield_process_id = 1;
+ initializeProcessBoundaryConditionsAndSourceTerms(
+ *_local_to_global_index_map_single_component, phasefield_process_id);
}
template <int DisplacementDim>
@@ -307,7 +306,7 @@ void PhaseFieldProcess<DisplacementDim>::assembleWithJacobianConcreteProcess(
else
{
// For the staggered scheme
- if (_coupled_solutions->process_id == 0)
+ if (_coupled_solutions->process_id == 1)
{
DBUG(
"Assemble the Jacobian equations of phase field in "
@@ -319,10 +318,11 @@ void PhaseFieldProcess<DisplacementDim>::assembleWithJacobianConcreteProcess(
"Assemble the Jacobian equations of deformation in "
"PhaseFieldProcess for the staggered scheme.");
}
- dof_tables.emplace_back(*_local_to_global_index_map_single_component);
dof_tables.emplace_back(*_local_to_global_index_map);
+ dof_tables.emplace_back(*_local_to_global_index_map_single_component);
}
// Call global assembler for each local assembly item.
+
GlobalExecutor::executeMemberDereferenced(
_global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
_local_assemblers, dof_tables, t, x, xdot, dxdot_dx, dx_dx, M, K, b,
@@ -359,16 +359,12 @@ template <int DisplacementDim>
void PhaseFieldProcess<DisplacementDim>::postNonLinearSolverConcreteProcess(
GlobalVector const& x, const double t, const int process_id)
{
- if (!_use_monolithic_scheme && process_id == 0)
- {
- return;
- }
-
DBUG("PostNonLinearSolver PhaseFieldProcess.");
// Calculate strain, stress or other internal variables of mechanics.
GlobalExecutor::executeMemberOnDereferenced(
&LocalAssemblerInterface::postNonLinearSolver, _local_assemblers,
getDOFTable(process_id), x, t, _use_monolithic_scheme);
}
+
} // namespace PhaseField
} // namespace ProcessLib
diff --git a/ProcessLib/PhaseField/PhaseFieldProcess.h b/ProcessLib/PhaseField/PhaseFieldProcess.h
index ef31089b2b5bd7b3e3182609c1299a6412d4f627..009cb18506cbb8b51d94625b6e1aed3aa00c25d3 100644
--- a/ProcessLib/PhaseField/PhaseFieldProcess.h
+++ b/ProcessLib/PhaseField/PhaseFieldProcess.h
@@ -18,7 +18,6 @@ namespace ProcessLib
{
namespace PhaseField
{
-
/**
* \brief A class to simulate mechanical fracturing process
* using phase-field approach in solids described by
@@ -87,25 +86,25 @@ private:
MeshLib::Mesh const& mesh,
unsigned const integration_order) override;
- void assembleConcreteProcess(
- const double t, GlobalVector const& x, GlobalMatrix& M, GlobalMatrix& K,
- GlobalVector& b) override;
+ void assembleConcreteProcess(const double t, GlobalVector const& x,
+ GlobalMatrix& M, GlobalMatrix& K,
+ GlobalVector& b) override;
void assembleWithJacobianConcreteProcess(
const double t, GlobalVector const& x, GlobalVector const& xdot,
const double dxdot_dx, const double dx_dx, GlobalMatrix& M,
GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) override;
- void preTimestepConcreteProcess(
- GlobalVector const& x, double const t, double const dt,
- const int process_id) override;
+ void preTimestepConcreteProcess(GlobalVector const& x, double const t,
+ double const dt,
+ const int process_id) override;
void postTimestepConcreteProcess(GlobalVector const& x,
int const process_id) override;
void postNonLinearSolverConcreteProcess(GlobalVector const& x,
const double t,
- int const processs_id) override;
+ int const process_id) override;
private:
PhaseFieldProcessData<DisplacementDim> _process_data;
@@ -116,8 +115,16 @@ private:
_local_to_global_index_map_single_component;
/// Sparsity pattern for the phase field equation, and it is initialized
- // only if the staggered scheme is used.
+ /// only if the staggered scheme is used.
GlobalSparsityPattern _sparsity_pattern_with_single_component;
+
+ /// Check whether the process represented by \c process_id is/has
+ /// mechanical process. In the present implementation, the mechanical
+ /// process has process_id == 0 in the staggered scheme.
+ bool hasMechanicalProcess(int const process_id) const
+ {
+ return _use_monolithic_scheme || process_id == 0;
+ }
};
extern template class PhaseFieldProcess<2>;