diff --git a/ProcessLib/PhaseField/PhaseFieldFEM-impl.h b/ProcessLib/PhaseField/PhaseFieldFEM-impl.h
index 25db542f1081d1207e429a378e5b0956925a095f..812a4b36b1e6c5d554529354cec95ef3e29a68aa 100644
--- a/ProcessLib/PhaseField/PhaseFieldFEM-impl.h
+++ b/ProcessLib/PhaseField/PhaseFieldFEM-impl.h
@@ -14,6 +14,187 @@ namespace ProcessLib
{
namespace PhaseField
{
+template <typename ShapeFunction, typename IntegrationMethod,
+ int DisplacementDim>
+void PhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
+ DisplacementDim>::
+ assembleWithJacobian(double const t, std::vector<double> const& local_x,
+ 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_rhs_data,
+ std::vector<double>& local_Jac_data)
+{
+ auto const local_matrix_size = local_x.size();
+ assert(local_matrix_size == phasefield_size + displacement_size);
+
+ auto d = Eigen::Map<
+ typename ShapeMatricesType::template VectorType<phasefield_size> const>(
+ local_x.data() + phasefield_index, phasefield_size);
+
+ auto u = Eigen::Map<typename ShapeMatricesType::template VectorType<
+ displacement_size> const>(local_x.data() + displacement_index,
+ displacement_size);
+
+ auto d_dot = Eigen::Map<
+ typename ShapeMatricesType::template VectorType<phasefield_size> const>(
+ local_xdot.data() + phasefield_index, phasefield_size);
+
+ auto local_Jac = MathLib::createZeroedMatrix<JacobianMatrix>(
+ local_Jac_data, local_matrix_size, local_matrix_size);
+
+ auto local_rhs = MathLib::createZeroedVector<RhsVector>(local_rhs_data,
+ local_matrix_size);
+
+ typename ShapeMatricesType::template MatrixType<displacement_size,
+ phasefield_size>
+ Kud;
+ Kud.setZero(displacement_size, phasefield_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());
+
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ x_position.setIntegrationPoint(ip);
+ 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 const& C_tensile = _ip_data[ip].C_tensile;
+ auto const& C_compressive = _ip_data[ip].C_compressive;
+
+ auto& eps = _ip_data[ip].eps;
+ auto const& strain_energy_tensile = _ip_data[ip].strain_energy_tensile;
+ auto const& sigma_tensile = _ip_data[ip].sigma_tensile;
+
+ auto& history_variable = _ip_data[ip].history_variable;
+ auto& history_variable_prev = _ip_data[ip].history_variable_prev;
+
+ auto const& sigma_real = _ip_data[ip].sigma_real;
+
+ // Kdd_1 defines one term which both used in Kdd and local_rhs for
+ // phase field
+ double const gc = _process_data.crack_resistance(t, x_position)[0];
+ double const ls = _process_data.crack_length_scale(t, x_position)[0];
+ typename ShapeMatricesType::NodalMatrixType const Kdd_1 =
+ dNdx.transpose() * 2 * gc * ls * dNdx;
+
+ //
+ // displacement equation, displacement part
+ //
+ 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;
+ eps.noalias() = B * u;
+ _ip_data[ip].updateConstitutiveRelation(t, x_position, dt, u,
+ degradation);
+
+ local_Jac
+ .template block<displacement_size, displacement_size>(
+ displacement_index, displacement_index)
+ .noalias() +=
+ B.transpose() * (degradation * C_tensile + C_compressive) * B * w;
+
+ 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;
+ local_rhs.template block<displacement_size, 1>(displacement_index, 0)
+ .noalias() -=
+ (B.transpose() * sigma_real - N_u.transpose() * rho_sr * b) * w;
+
+ //
+ // displacement equation, phasefield part
+ //
+
+ // Temporary storage used in the Kud and for potential reuse in the
+ // Kdu matrix.
+ auto const Kud_ip_contribution =
+ (B.transpose() * 2 * d_ip * sigma_tensile * N * w).eval();
+
+ Kud.noalias() += Kud_ip_contribution;
+
+ if (history_variable_prev < strain_energy_tensile)
+ {
+ history_variable = strain_energy_tensile;
+ Kdu.noalias() += Kud_ip_contribution.transpose();
+ }
+ else
+ {
+ history_variable = history_variable_prev;
+ }
+
+ //
+ // phasefield equation, phasefield part.
+ //
+ Kdd.noalias() += (Kdd_1 + N.transpose() * 2 * history_variable * N +
+ N.transpose() * 0.5 * gc / ls * N) *
+ w;
+ double const M = _process_data.kinetic_coefficient(t, x_position)[0];
+ double const d_dot_ip = N.dot(d_dot);
+
+ local_rhs.template segment<phasefield_size>(phasefield_index)
+ .noalias() -= (N.transpose() * d_dot_ip / M + Kdd_1 * d +
+ N.transpose() * d_ip * 2 * history_variable -
+ N.transpose() * 0.5 * gc / ls * (1 - d_ip)) *
+ w;
+
+ Ddd.noalias() += N.transpose() / M * N * w;
+ }
+ // displacement equation, phasefield part
+ local_Jac
+ .template block<displacement_size, phasefield_size>(displacement_index,
+ phasefield_index)
+ .noalias() += Kud;
+
+ // phasefield equation, phasefield part.
+ local_Jac
+ .template block<phasefield_size, phasefield_size>(phasefield_index,
+ phasefield_index)
+ .noalias() += Kdd + Ddd / dt;
+
+ // phasefield equation, displacement part.
+ local_Jac
+ .template block<phasefield_size, displacement_size>(phasefield_index,
+ displacement_index)
+ .noalias() += Kdu;
+}
+
template <typename ShapeFunction, typename IntegrationMethod,
int DisplacementDim>
void PhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
diff --git a/ProcessLib/PhaseField/PhaseFieldFEM.h b/ProcessLib/PhaseField/PhaseFieldFEM.h
index 976fa94a78bdc2f04f913e00efd6576883f6dbaa..7486f72abf1b544f0d71b74995ee56cd3712501d 100644
--- a/ProcessLib/PhaseField/PhaseFieldFEM.h
+++ b/ProcessLib/PhaseField/PhaseFieldFEM.h
@@ -190,180 +190,7 @@ public:
std::vector<double>& /*local_M_data*/,
std::vector<double>& /*local_K_data*/,
std::vector<double>& local_rhs_data,
- std::vector<double>& local_Jac_data) override
- {
- auto const local_matrix_size = local_x.size();
- assert(local_matrix_size == phasefield_size + displacement_size);
-
- auto d = Eigen::Map<typename ShapeMatricesType::template VectorType<
- phasefield_size> const>(local_x.data() + phasefield_index,
- phasefield_size);
-
- auto u = Eigen::Map<typename ShapeMatricesType::template VectorType<
- displacement_size> const>(local_x.data() + displacement_index,
- displacement_size);
-
- auto d_dot = Eigen::Map<typename ShapeMatricesType::template VectorType<
- phasefield_size> const>(local_xdot.data() + phasefield_index,
- phasefield_size);
-
- auto local_Jac = MathLib::createZeroedMatrix<JacobianMatrix>(
- local_Jac_data, local_matrix_size, local_matrix_size);
-
- auto local_rhs = MathLib::createZeroedVector<RhsVector>(
- local_rhs_data, local_matrix_size);
-
- typename ShapeMatricesType::template MatrixType<displacement_size,
- phasefield_size>
- Kud;
- Kud.setZero(displacement_size, phasefield_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());
-
- unsigned const n_integration_points =
- _integration_method.getNumberOfPoints();
- for (unsigned ip = 0; ip < n_integration_points; ip++)
- {
- x_position.setIntegrationPoint(ip);
- 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 const& C_tensile = _ip_data[ip].C_tensile;
- auto const& C_compressive = _ip_data[ip].C_compressive;
-
- auto& eps = _ip_data[ip].eps;
- auto const& strain_energy_tensile =
- _ip_data[ip].strain_energy_tensile;
- auto const& sigma_tensile = _ip_data[ip].sigma_tensile;
-
- auto& history_variable = _ip_data[ip].history_variable;
- auto& history_variable_prev = _ip_data[ip].history_variable_prev;
-
- auto const& sigma_real = _ip_data[ip].sigma_real;
-
- // Kdd_1 defines one term which both used in Kdd and local_rhs for
- // phase field
- double const gc = _process_data.crack_resistance(t, x_position)[0];
- double const ls =
- _process_data.crack_length_scale(t, x_position)[0];
- typename ShapeMatricesType::NodalMatrixType const Kdd_1 =
- dNdx.transpose() * 2 * gc * ls * dNdx;
-
- //
- // displacement equation, displacement part
- //
- 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;
- eps.noalias() = B * u;
- _ip_data[ip].updateConstitutiveRelation(t, x_position, dt, u,
- degradation);
-
- local_Jac
- .template block<displacement_size, displacement_size>(
- displacement_index, displacement_index)
- .noalias() += B.transpose() *
- (degradation * C_tensile + C_compressive) * B * w;
-
- 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;
- local_rhs
- .template block<displacement_size, 1>(displacement_index, 0)
- .noalias() -=
- (B.transpose() * sigma_real - N_u.transpose() * rho_sr * b) * w;
-
- //
- // displacement equation, phasefield part
- //
-
- // Temporary storage used in the Kud and for potential reuse in the
- // Kdu matrix.
- auto const Kud_ip_contribution =
- (B.transpose() * 2 * d_ip * sigma_tensile * N * w).eval();
-
- Kud.noalias() += Kud_ip_contribution;
-
- if (history_variable_prev < strain_energy_tensile)
- {
- history_variable = strain_energy_tensile;
- Kdu.noalias() += Kud_ip_contribution.transpose();
- }
- else
- {
- history_variable = history_variable_prev;
- }
-
- //
- // phasefield equation, phasefield part.
- //
- Kdd.noalias() += (Kdd_1 + N.transpose() * 2 * history_variable * N +
- N.transpose() * 0.5 * gc / ls * N) *
- w;
- double const M =
- _process_data.kinetic_coefficient(t, x_position)[0];
- double const d_dot_ip = N.dot(d_dot);
-
- local_rhs.template segment<phasefield_size>(phasefield_index)
- .noalias() -= (N.transpose() * d_dot_ip / M + Kdd_1 * d +
- N.transpose() * d_ip * 2 * history_variable -
- N.transpose() * 0.5 * gc / ls * (1 - d_ip)) *
- w;
-
- Ddd.noalias() += N.transpose() / M * N * w;
- }
- // displacement equation, phasefield part
- local_Jac
- .template block<displacement_size, phasefield_size>(
- displacement_index, phasefield_index)
- .noalias() += Kud;
-
- // phasefield equation, phasefield part.
- local_Jac
- .template block<phasefield_size, phasefield_size>(phasefield_index,
- phasefield_index)
- .noalias() += Kdd + Ddd / dt;
-
- // phasefield equation, displacement part.
- local_Jac
- .template block<phasefield_size, displacement_size>(
- phasefield_index, displacement_index)
- .noalias() += Kdu;
- }
+ std::vector<double>& local_Jac_data) override;
void assembleWithJacobianForStaggeredScheme(
double const t, std::vector<double> const& local_xdot,