OGS
HTProcess.cpp
Go to the documentation of this file.
1 
11 #include "HTProcess.h"
12 
13 #include <cassert>
14 
15 #include "MonolithicHTFEM.h"
20 #include "StaggeredHTFEM.h"
21 
22 namespace ProcessLib
23 {
24 namespace HT
25 {
27  std::string name,
28  MeshLib::Mesh& mesh,
29  std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
30  std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
31  unsigned const integration_order,
32  std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
33  process_variables,
34  HTProcessData&& process_data,
35  SecondaryVariableCollection&& secondary_variables,
36  bool const use_monolithic_scheme,
37  std::unique_ptr<ProcessLib::SurfaceFluxData>&& surfaceflux)
38  : Process(std::move(name), mesh, std::move(jacobian_assembler), parameters,
39  integration_order, std::move(process_variables),
40  std::move(secondary_variables), use_monolithic_scheme),
41  _process_data(std::move(process_data)),
42  _surfaceflux(std::move(surfaceflux))
43 {
44 }
45 
47  NumLib::LocalToGlobalIndexMap const& dof_table,
48  MeshLib::Mesh const& mesh,
49  unsigned const integration_order)
50 {
51  // For the staggered scheme, both processes are assumed to use the same
52  // element order. Therefore the order of shape function can be fetched from
53  // any set of the sets of process variables of the coupled processes. Here,
54  // we take the one from the first process by setting process_id = 0.
55  const int process_id = 0;
56  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
57 
59  {
60  ProcessLib::createLocalAssemblers<MonolithicHTFEM>(
61  mesh.getDimension(), mesh.getElements(), dof_table,
63  mesh.isAxiallySymmetric(), integration_order, _process_data);
64  }
65  else
66  {
67  ProcessLib::createLocalAssemblers<StaggeredHTFEM>(
68  mesh.getDimension(), mesh.getElements(), dof_table,
70  mesh.isAxiallySymmetric(), integration_order, _process_data);
71  }
72 
74  "darcy_velocity",
78 }
79 
80 void HTProcess::assembleConcreteProcess(const double t, double const dt,
81  std::vector<GlobalVector*> const& x,
82  std::vector<GlobalVector*> const& xdot,
83  int const process_id, GlobalMatrix& M,
85 {
86  std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
87  dof_tables;
89  {
90  DBUG("Assemble HTProcess.");
91  dof_tables.emplace_back(*_local_to_global_index_map);
92  }
93  else
94  {
95  if (process_id == _process_data.heat_transport_process_id)
96  {
97  DBUG(
98  "Assemble the equations of heat transport process within "
99  "HTProcess.");
100  }
101  else
102  {
103  DBUG(
104  "Assemble the equations of single phase fully saturated "
105  "fluid flow process within HTProcess.");
106  }
107  dof_tables.emplace_back(*_local_to_global_index_map);
108  dof_tables.emplace_back(*_local_to_global_index_map);
109  }
110 
111  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
112  // Call global assembler for each local assembly item.
115  pv.getActiveElementIDs(), dof_tables, t, dt, x, xdot, process_id, M, K,
116  b);
117 }
118 
120  const double t, double const dt, std::vector<GlobalVector*> const& x,
121  std::vector<GlobalVector*> const& xdot, const double dxdot_dx,
122  const double dx_dx, int const process_id, GlobalMatrix& M, GlobalMatrix& K,
123  GlobalVector& b, GlobalMatrix& Jac)
124 {
125  DBUG("AssembleWithJacobian HTProcess.");
126 
127  std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
128  dof_tables;
130  {
131  dof_tables.emplace_back(std::ref(*_local_to_global_index_map));
132  }
133  else
134  {
135  dof_tables.emplace_back(std::ref(*_local_to_global_index_map));
136  dof_tables.emplace_back(std::ref(*_local_to_global_index_map));
137  }
138 
139  // Call global assembler for each local assembly item.
140  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
143  _local_assemblers, pv.getActiveElementIDs(), dof_tables, t, dt, x, xdot,
144  dxdot_dx, dx_dx, process_id, M, K, b, Jac);
145 }
146 
148  int const process_id)
149 {
150  DBUG("Set the coupled term for the staggered scheme to local assembers.");
151 
152  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
156 }
157 
158 std::tuple<NumLib::LocalToGlobalIndexMap*, bool>
160 {
162  {
163  // For single-variable-single-component processes reuse the existing DOF
164  // table.
165  const bool manage_storage = false;
166  return std::make_tuple(_local_to_global_index_map.get(),
167  manage_storage);
168  }
169 
170  // Otherwise construct a new DOF table.
171  std::vector<MeshLib::MeshSubset> all_mesh_subsets_single_component{
173 
174  const bool manage_storage = true;
175  return std::make_tuple(new NumLib::LocalToGlobalIndexMap(
176  std::move(all_mesh_subsets_single_component),
177  // by location order is needed for output
179  manage_storage);
180 }
181 
182 Eigen::Vector3d HTProcess::getFlux(std::size_t element_id,
183  MathLib::Point3d const& p,
184  double const t,
185  std::vector<GlobalVector*> const& x) const
186 {
187  // fetch local_x from primary variable
188  std::vector<GlobalIndexType> indices_cache;
189  auto const r_c_indices = NumLib::getRowColumnIndices(
190  element_id, *_local_to_global_index_map, indices_cache);
191  std::vector<std::vector<GlobalIndexType>> indices_of_all_coupled_processes{
192  x.size(), r_c_indices.rows};
193  auto const local_x =
194  getCoupledLocalSolutions(x, indices_of_all_coupled_processes);
195 
196  return _local_assemblers[element_id]->getFlux(p, t, local_x);
197 }
198 
199 // this is almost a copy of the implementation in the GroundwaterFlow
200 void HTProcess::postTimestepConcreteProcess(std::vector<GlobalVector*> const& x,
201  const double t,
202  const double /*delta_t*/,
203  int const process_id)
204 {
205  // For the monolithic scheme, process_id is always zero.
206  if (_use_monolithic_scheme && process_id != 0)
207  {
208  OGS_FATAL(
209  "The condition of process_id = 0 must be satisfied for monolithic "
210  "HTProcess, which is a single process.");
211  }
212  if (!_use_monolithic_scheme &&
213  process_id != _process_data.hydraulic_process_id)
214  {
215  DBUG("This is the thermal part of the staggered HTProcess.");
216  return;
217  }
218  if (!_surfaceflux) // computing the surfaceflux is optional
219  {
220  return;
221  }
222 
223  ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
224 
225  _surfaceflux->integrate(x, t, *this, process_id, _integration_order, _mesh,
226  pv.getActiveElementIDs());
227 }
228 } // namespace HT
229 } // namespace ProcessLib
#define OGS_FATAL(...)
Definition: Error.h:26
void DBUG(char const *fmt, Args const &... args)
Definition: Logging.h:27
Global vector based on Eigen vector.
Definition: EigenVector.h:26
bool isAxiallySymmetric() const
Definition: Mesh.h:126
unsigned getDimension() const
Returns the dimension of the mesh (determined by the maximum dimension over all elements).
Definition: Mesh.h:71
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition: Mesh.h:98
virtual std::vector< double > const & getIntPtDarcyVelocity(const double, std::vector< GlobalVector * > const &x, std::vector< NumLib::LocalToGlobalIndexMap const * > const &, std::vector< double > &) const =0
void setStaggeredCoupledSolutions(std::size_t const, CoupledSolutionsForStaggeredScheme *const coupling_term)
void postTimestepConcreteProcess(std::vector< GlobalVector * > const &x, const double t, const double delta_t, int const process_id) override
Definition: HTProcess.cpp:200
void setCoupledTermForTheStaggeredSchemeToLocalAssemblers(int const process_id) override
Definition: HTProcess.cpp:147
void initializeConcreteProcess(NumLib::LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh, unsigned const integration_order) override
Process specific initialization called by initialize().
Definition: HTProcess.cpp:46
void assembleConcreteProcess(const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &xdot, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b) override
Definition: HTProcess.cpp:80
std::tuple< NumLib::LocalToGlobalIndexMap *, bool > getDOFTableForExtrapolatorData() const override
Definition: HTProcess.cpp:159
std::unique_ptr< ProcessLib::SurfaceFluxData > _surfaceflux
Definition: HTProcess.h:117
Eigen::Vector3d getFlux(std::size_t element_id, MathLib::Point3d const &p, double const t, std::vector< GlobalVector * > const &x) const override
Definition: HTProcess.cpp:182
HTProcessData _process_data
Definition: HTProcess.h:113
HTProcess(std::string name, MeshLib::Mesh &mesh, std::unique_ptr< ProcessLib::AbstractJacobianAssembler > &&jacobian_assembler, std::vector< std::unique_ptr< ParameterLib::ParameterBase >> const &parameters, unsigned const integration_order, std::vector< std::vector< std::reference_wrapper< ProcessVariable >>> &&process_variables, HTProcessData &&process_data, SecondaryVariableCollection &&secondary_variables, bool const use_monolithic_scheme, std::unique_ptr< ProcessLib::SurfaceFluxData > &&surfaceflux)
Definition: HTProcess.cpp:26
void assembleWithJacobianConcreteProcess(const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &xdot, const double dxdot_dx, const double dx_dx, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, GlobalMatrix &Jac) override
Definition: HTProcess.cpp:119
std::vector< std::unique_ptr< HTLocalAssemblerInterface > > _local_assemblers
Definition: HTProcess.h:115
unsigned getShapeFunctionOrder() const
std::vector< std::size_t > const & getActiveElementIDs() const
NumLib::Extrapolator & getExtrapolator() const
Definition: Process.h:185
std::unique_ptr< MeshLib::MeshSubset const > _mesh_subset_all_nodes
Definition: Process.h:327
MeshLib::Mesh & _mesh
Definition: Process.h:326
SecondaryVariableCollection _secondary_variables
Definition: Process.h:331
CoupledSolutionsForStaggeredScheme * _coupled_solutions
Definition: Process.h:339
std::vector< std::reference_wrapper< ProcessVariable > > const & getProcessVariables(const int process_id) const
Definition: Process.h:145
VectorMatrixAssembler _global_assembler
Definition: Process.h:333
unsigned const _integration_order
Definition: Process.h:344
std::unique_ptr< NumLib::LocalToGlobalIndexMap > _local_to_global_index_map
Definition: Process.h:329
const bool _use_monolithic_scheme
Definition: Process.h:335
Handles configuration of several secondary variables from the project file.
void addSecondaryVariable(std::string const &internal_name, SecondaryVariableFunctions &&fcts)
void assembleWithJacobian(std::size_t const mesh_item_id, LocalAssemblerInterface &local_assembler, std::vector< std::reference_wrapper< NumLib::LocalToGlobalIndexMap >> const &dof_tables, const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &xdot, const double dxdot_dx, const double dx_dx, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b, GlobalMatrix &Jac)
void assemble(std::size_t const mesh_item_id, LocalAssemblerInterface &local_assembler, std::vector< std::reference_wrapper< NumLib::LocalToGlobalIndexMap >> const &dof_tables, double const t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &xdot, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b)
NumLib::LocalToGlobalIndexMap::RowColumnIndices getRowColumnIndices(std::size_t const id, NumLib::LocalToGlobalIndexMap const &dof_table, std::vector< GlobalIndexType > &indices)
@ BY_LOCATION
Ordering data by spatial location.
std::vector< double > getCoupledLocalSolutions(std::vector< GlobalVector * > const &global_solutions, std::vector< std::vector< GlobalIndexType >> const &indices)
SecondaryVariableFunctions makeExtrapolator(const unsigned num_components, NumLib::Extrapolator &extrapolator, LocalAssemblerCollection const &local_assemblers, typename NumLib::ExtrapolatableLocalAssemblerCollection< LocalAssemblerCollection >::IntegrationPointValuesMethod integration_point_values_method)
static void executeSelectedMemberOnDereferenced(Method method, Container const &container, std::vector< std::size_t > const &active_container_ids, Args &&... args)
static void executeSelectedMemberDereferenced(Object &object, Method method, Container const &container, std::vector< std::size_t > const &active_container_ids, Args &&... args)