OGS
HydroMechanicsLocalAssemblerMatrixNearFracture-impl.h
Go to the documentation of this file.
1
11#pragma once
12
15
16namespace ProcessLib
17{
18namespace LIE
19{
20namespace HydroMechanics
21{
22template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
23 int DisplacementDim>
24HydroMechanicsLocalAssemblerMatrixNearFracture<
25 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
26 HydroMechanicsLocalAssemblerMatrixNearFracture(
27 MeshLib::Element const& e,
28 std::size_t const n_variables,
29 std::size_t const local_matrix_size,
30 std::vector<unsigned> const& dofIndex_to_localIndex,
31 NumLib::GenericIntegrationMethod const& integration_method,
32 bool const is_axially_symmetric,
35 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>(
36 e, n_variables, local_matrix_size, dofIndex_to_localIndex,
37 integration_method, is_axially_symmetric, process_data),
38 _e_center_coords(getCenterOfGravity(e).data())
39{
40 for (auto fid : process_data.vec_ele_connected_fractureIDs[e.getID()])
41 {
42 _fracID_to_local.insert({fid, _fracture_props.size()});
43 _fracture_props.push_back(&_process_data.fracture_properties[fid]);
44 }
45
46 // TODO(naumov) Junctions not yet implemented.
47}
48
49template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
50 int DisplacementDim>
51void HydroMechanicsLocalAssemblerMatrixNearFracture<
52 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
53 assembleWithJacobianConcrete(double const t, double const dt,
54 Eigen::VectorXd const& local_x,
55 Eigen::VectorXd const& local_x_prev,
56 Eigen::VectorXd& local_b,
57 Eigen::MatrixXd& local_J)
58{
59 auto p = const_cast<Eigen::VectorXd&>(local_x).segment(pressure_index,
60 pressure_size);
61 auto p_prev = const_cast<Eigen::VectorXd&>(local_x_prev)
62 .segment(pressure_index, pressure_size);
63 if (_process_data.deactivate_matrix_in_flow)
64 {
65 Base::setPressureOfInactiveNodes(t, p);
66 }
67 auto const u = local_x.segment(displacement_index, displacement_size);
68 auto const u_prev =
69 local_x_prev.segment(displacement_index, displacement_size);
70
71 auto rhs_p = local_b.segment(pressure_index, pressure_size);
72 auto rhs_u = local_b.segment(displacement_index, displacement_size);
73
74 auto J_pp = local_J.block(pressure_index, pressure_index, pressure_size,
75 pressure_size);
76 auto J_pu = local_J.block(pressure_index, displacement_index, pressure_size,
77 displacement_size);
78 auto J_up = local_J.block(displacement_index, pressure_index,
79 displacement_size, pressure_size);
80 auto J_uu = local_J.block(displacement_index, displacement_index,
81 displacement_size, displacement_size);
82
83 // levelset value of the element
84 // remark: this assumes the levelset function is uniform within an element
85 std::vector<double> levelsets = uGlobalEnrichments(
86 _fracture_props, _junction_props, _fracID_to_local, _e_center_coords);
87 double const ele_levelset = levelsets[0]; // single fracture
88
89 if (ele_levelset == 0)
90 {
91 // no DoF exists for displacement jumps. do the normal assembly
92 Base::assembleBlockMatricesWithJacobian(
93 t, dt, p, p_prev, u, u_prev, rhs_p, rhs_u, J_pp, J_pu, J_uu, J_up);
94 return;
95 }
96
97 // Displacement jumps should be taken into account
98
99 // compute true displacements
100 auto const g = local_x.segment(displacement_jump_index, displacement_size);
101 auto const g_prev =
102 local_x_prev.segment(displacement_jump_index, displacement_size);
103 Eigen::VectorXd const total_u = u + ele_levelset * g;
104 Eigen::VectorXd const total_u_prev = u_prev + ele_levelset * g_prev;
105
106 // evaluate residuals and Jacobians for pressure and displacements
107 Base::assembleBlockMatricesWithJacobian(t, dt, p, p_prev, total_u,
108 total_u_prev, rhs_p, rhs_u, J_pp,
109 J_pu, J_uu, J_up);
110
111 // compute residuals and Jacobians for displacement jumps
112 auto rhs_g = local_b.segment(displacement_jump_index, displacement_size);
113 auto J_pg = local_J.block(pressure_index, displacement_jump_index,
114 pressure_size, displacement_size);
115 auto J_ug = local_J.block(displacement_index, displacement_jump_index,
116 displacement_size, displacement_size);
117 auto J_gp = local_J.block(displacement_jump_index, pressure_index,
118 displacement_size, pressure_size);
119 auto J_gu = local_J.block(displacement_jump_index, displacement_index,
120 displacement_size, displacement_size);
121 auto J_gg = local_J.block(displacement_jump_index, displacement_jump_index,
122 displacement_size, displacement_size);
123
124 rhs_g = ele_levelset * rhs_u;
125 J_pg = ele_levelset * J_pu;
126 J_ug = ele_levelset * J_uu;
127 J_gp = ele_levelset * J_up;
128 J_gu = ele_levelset * J_uu;
129 J_gg = ele_levelset * ele_levelset * J_uu;
130}
131
132template <typename ShapeFunctionDisplacement, typename ShapeFunctionPressure,
133 int DisplacementDim>
135 ShapeFunctionDisplacement, ShapeFunctionPressure, DisplacementDim>::
136 postTimestepConcreteWithVector(double const t, double const dt,
137 Eigen::VectorXd const& local_x)
138{
139 auto p = const_cast<Eigen::VectorXd&>(local_x).segment(pressure_index,
140 pressure_size);
141 if (_process_data.deactivate_matrix_in_flow)
142 {
143 Base::setPressureOfInactiveNodes(t, p);
144 }
145 auto u = local_x.segment(displacement_index, displacement_size);
146
147 // levelset value of the element
148 // remark: this assumes the levelset function is uniform within an element
149 std::vector<double> levelsets = uGlobalEnrichments(
150 _fracture_props, _junction_props, _fracID_to_local, _e_center_coords);
151 double const ele_levelset = levelsets[0]; // single fracture
152
153 if (ele_levelset == 0)
154 {
155 // no DoF exists for displacement jumps. do the normal assembly
156 Base::postTimestepConcreteWithBlockVectors(t, dt, p, u);
157 return;
158 }
159
160 // Displacement jumps should be taken into account
161
162 // compute true displacements
163 auto const g = local_x.segment(displacement_jump_index, displacement_size);
164 Eigen::VectorXd const total_u = u + ele_levelset * g;
165
166 // evaluate residuals and Jacobians for pressure and displacements
167 Base::postTimestepConcreteWithBlockVectors(t, dt, p, total_u);
168}
169
170} // namespace HydroMechanics
171} // namespace LIE
172} // namespace ProcessLib
std::size_t getID() const
Returns the ID of the element.
Definition Element.h:91
std::vector< double > uGlobalEnrichments(std::vector< FractureProperty * > const &frac_props, std::vector< JunctionProperty * > const &junction_props, std::unordered_map< int, int > const &fracID_to_local, Eigen::Vector3d const &x)