PolyDiM
C++ library for POLYtopal DIscretization Methods
Loading...
Searching...
No Matches
Assembler_Utilities.hpp
Go to the documentation of this file.
1// _LICENSE_HEADER_
2//
3// Copyright (C) 2019 - 2025.
4// Terms register on the GPL-3.0 license.
5//
6// This file can be redistributed and/or modified under the license terms.
7//
8// See top level LICENSE file for more details.
9//
10// This file can be used citing references in CITATION.cff file.
11
12#ifndef __PDETOOLS_ASSEMBLER_AssemblerUtilities_HPP
13#define __PDETOOLS_ASSEMBLER_AssemblerUtilities_HPP
14
15#include "DOFsManager.hpp"
16#include <vector>
17
18namespace Polydim
19{
20namespace PDETools
21{
22namespace Assembler_Utilities
23{
25{
26 std::vector<std::reference_wrapper<const Polydim::PDETools::DOFs::DOFsManager::DOFsData>> dofs_data;
27 std::vector<size_t> local_offsets;
28 std::vector<size_t> global_offsets_DOFs;
29 std::vector<size_t> global_offsets_Strongs;
30};
31
32struct count_dofs_data final
33{
35 unsigned int num_total_dofs;
36 unsigned int num_total_strong;
37 std::vector<size_t> offsets_DOFs;
38 std::vector<size_t> offsets_Strongs;
39};
40
42{
43 unsigned int num_total_dofs;
44 std::vector<size_t> offsets_DOFs;
45};
46
47inline Polydim::PDETools::Assembler_Utilities::count_dofs_data count_dofs(const std::vector<Polydim::PDETools::DOFs::DOFsManager::DOFsData> &dofs_data)
48{
50 data.num_total_dofs = dofs_data[0].NumberDOFs;
51 data.num_total_strong = dofs_data[0].NumberStrongs;
52 data.num_total_boundary_dofs = dofs_data[0].NumberBoundaryDOFs;
53
54 const unsigned int numDOFHandler = dofs_data.size();
55 data.offsets_DOFs.resize(numDOFHandler);
56 data.offsets_Strongs.resize(numDOFHandler);
57 data.offsets_DOFs[0] = 0;
58 data.offsets_Strongs[0] = 0;
59
60 for (unsigned int i = 1; i < numDOFHandler; i++)
61 {
62 data.num_total_dofs += dofs_data[i].NumberDOFs;
63 data.num_total_strong += dofs_data[i].NumberStrongs;
64 data.num_total_boundary_dofs += dofs_data[i].NumberBoundaryDOFs;
65
66 data.offsets_DOFs[i] = data.offsets_DOFs[i - 1] + dofs_data[i - 1].NumberDOFs;
67 data.offsets_Strongs[i] = data.offsets_Strongs[i - 1] + dofs_data[i - 1].NumberStrongs;
68 }
69
70 return data;
71}
72// ***************************************************************************
73template <unsigned int dimension>
75 const unsigned int cell_index,
76 const std::vector<std::reference_wrapper<const Polydim::PDETools::DOFs::DOFsManager::DOFsData>> &dofs_data)
77{
79 data.num_total_dofs = dofs_data[0].get().CellsGlobalDOFs[dimension].at(cell_index).size();
80
81 const unsigned int numDOFHandler = dofs_data.size();
82 data.offsets_DOFs.resize(numDOFHandler);
83 data.offsets_DOFs[0] = 0;
84
85 for (unsigned int i = 1; i < numDOFHandler; i++)
86 {
87 data.num_total_dofs += dofs_data[i].get().CellsGlobalDOFs[dimension].at(cell_index).size();
88 data.offsets_DOFs[i] =
89 data.offsets_DOFs[i - 1] + dofs_data[i - 1].get().CellsGlobalDOFs[dimension].at(cell_index).size();
90 }
91
92 return data;
93}
94// ***************************************************************************
95template <unsigned int dimension>
97local_count_dofs(const unsigned int cell_index, const Polydim::PDETools::DOFs::DOFsManager::DOFsData &dofs_data)
98{
99 std::vector<std::reference_wrapper<const Polydim::PDETools::DOFs::DOFsManager::DOFsData>> dofs_data_ref;
100 dofs_data_ref.reserve(1);
101 dofs_data_ref.push_back(std::cref(static_cast<const Polydim::PDETools::DOFs::DOFsManager::DOFsData &>(dofs_data)));
102
103 return local_count_dofs<dimension>(cell_index, dofs_data_ref);
104}
105// ***************************************************************************
106template <unsigned int dimension>
108local_count_dofs(const unsigned int cell_index, const std::vector<Polydim::PDETools::DOFs::DOFsManager::DOFsData> &dofs_data)
109{
110 std::vector<std::reference_wrapper<const Polydim::PDETools::DOFs::DOFsManager::DOFsData>> dofs_data_ref;
111 dofs_data_ref.reserve(dofs_data.size());
112 for (unsigned int c = 0; c < dofs_data.size(); c++)
113 dofs_data_ref.push_back(std::cref(static_cast<const Polydim::PDETools::DOFs::DOFsManager::DOFsData &>(dofs_data.at(c))));
114
115 return local_count_dofs<dimension>(cell_index, dofs_data_ref);
116}
117// ***************************************************************************
118template <unsigned int dimension, typename global_solution_type>
120 const unsigned int cell_index,
121 const std::vector<std::reference_wrapper<const Polydim::PDETools::DOFs::DOFsManager::DOFsData>> &dofs_data,
122 const unsigned int &num_local_DOFs,
123 const std::vector<size_t> &local_offsets_DOFs,
124 const std::vector<size_t> &global_offsets_DOFs,
125 const std::vector<size_t> &global_offsets_Strongs,
126 const global_solution_type &global_solution_DOFs,
127 const global_solution_type &global_solution_Strongs)
128{
129 Eigen::VectorXd local_solution_dofs = Eigen::VectorXd::Zero(num_local_DOFs);
130
131 const unsigned int numDOFHandler = dofs_data.size();
132 for (unsigned int h = 0; h < numDOFHandler; h++)
133 {
134 const auto &dofs = dofs_data[h].get();
135 const auto &global_dof = dofs.CellsGlobalDOFs[dimension].at(cell_index);
136 for (unsigned int loc_i = 0; loc_i < global_dof.size(); ++loc_i)
137 {
138 const auto &global_dof_i = global_dof.at(loc_i);
139 const auto &local_dof_i =
140 dofs.CellsDOFs.at(global_dof_i.Dimension).at(global_dof_i.CellIndex).at(global_dof_i.DOFIndex);
141
142 switch (local_dof_i.Type)
143 {
145 local_solution_dofs[loc_i + local_offsets_DOFs[h]] =
146 global_solution_Strongs.GetValue(local_dof_i.Global_Index + global_offsets_Strongs[h]);
147 break;
149 local_solution_dofs[loc_i + local_offsets_DOFs[h]] =
150 global_solution_DOFs.GetValue(local_dof_i.Global_Index + global_offsets_DOFs[h]);
151 break;
152 default:
153 throw std::runtime_error("Unknown DOF Type");
154 }
155 }
156 }
157
158 return local_solution_dofs;
159}
160// ***************************************************************************
161template <unsigned int dimension, typename global_solution_type>
162inline Eigen::VectorXd global_solution_to_local_solution(const unsigned int cell_index,
164 const unsigned int &num_local_DOFs,
165 const std::vector<size_t> &local_offsets_DOFs,
166 const std::vector<size_t> &global_offsets_DOFs,
167 const std::vector<size_t> &global_offsets_Strongs,
168 const global_solution_type &global_solution_DOFs,
169 const global_solution_type &global_solution_Strongs)
170{
171 std::vector<std::reference_wrapper<const Polydim::PDETools::DOFs::DOFsManager::DOFsData>> dofs_data_ref;
172 dofs_data_ref.reserve(1);
173 dofs_data_ref.push_back(std::cref(static_cast<const Polydim::PDETools::DOFs::DOFsManager::DOFsData &>(dofs_data)));
174
175 return global_solution_to_local_solution<dimension, global_solution_type>(cell_index,
176 dofs_data_ref,
177 num_local_DOFs,
178 local_offsets_DOFs,
179 global_offsets_DOFs,
180 global_offsets_Strongs,
181 global_solution_DOFs,
182 global_solution_Strongs);
183}
184// ***************************************************************************
185template <unsigned int dimension, typename global_solution_type>
186inline Eigen::VectorXd global_solution_to_local_solution(const unsigned int cell_index,
187 const std::vector<Polydim::PDETools::DOFs::DOFsManager::DOFsData> &dofs_data,
188 const unsigned int &num_local_DOFs,
189 const std::vector<size_t> &local_offsets_DOFs,
190 const std::vector<size_t> &global_offsets_DOFs,
191 const std::vector<size_t> &global_offsets_Strongs,
192 const global_solution_type &global_solution_DOFs,
193 const global_solution_type &global_solution_Strongs)
194{
195 std::vector<std::reference_wrapper<const Polydim::PDETools::DOFs::DOFsManager::DOFsData>> dofs_data_ref;
196 dofs_data_ref.reserve(dofs_data.size());
197 for (unsigned int c = 0; c < dofs_data.size(); c++)
198 dofs_data_ref.push_back(std::cref(static_cast<const Polydim::PDETools::DOFs::DOFsManager::DOFsData &>(dofs_data.at(c))));
199
200 return global_solution_to_local_solution<dimension, global_solution_type>(cell_index,
201 dofs_data_ref,
202 num_local_DOFs,
203 local_offsets_DOFs,
204 global_offsets_DOFs,
205 global_offsets_Strongs,
206 global_solution_DOFs,
207 global_solution_Strongs);
208}
209// ***************************************************************************
210template <unsigned int dimension, typename local_lhs_type, typename local_rhs_type, typename global_lhs_type, typename global_rhs_type>
211void assemble_local_matrix_to_global_matrix(const unsigned int cell_index,
212 const local_matrix_to_global_matrix_dofs_data &test_functions_dofs_data,
213 const local_matrix_to_global_matrix_dofs_data &trial_functions_dofs_data,
214 const local_lhs_type &local_lhs,
215 const local_rhs_type &local_rhs,
216 global_lhs_type &global_lhs_DOFs,
217 global_lhs_type &global_lhs_Strongs,
218 global_rhs_type &global_rhs)
219{
220 for (size_t test_f = 0; test_f < test_functions_dofs_data.dofs_data.size(); ++test_f)
221 {
222 const auto &test_dofs_data = test_functions_dofs_data.dofs_data.at(test_f).get();
223 const auto &test_global_dofs = test_dofs_data.CellsGlobalDOFs.at(dimension).at(cell_index);
224 const auto test_global_offset_DOFs = test_functions_dofs_data.global_offsets_DOFs[test_f];
225 const auto test_local_offset = test_functions_dofs_data.local_offsets[test_f];
226
227 for (size_t test_loc_i = 0; test_loc_i < test_global_dofs.size(); ++test_loc_i)
228 {
229 const auto global_dof_i = test_global_dofs.at(test_loc_i);
230 const auto local_dof_i =
231 test_dofs_data.CellsDOFs.at(global_dof_i.Dimension).at(global_dof_i.CellIndex).at(global_dof_i.DOFIndex);
232
233 switch (local_dof_i.Type)
234 {
236 continue;
238 break;
239 default:
240 throw std::runtime_error("Unknown DOF Type");
241 }
242
243 const unsigned int global_index_i = local_dof_i.Global_Index + test_global_offset_DOFs;
244
245 global_rhs.AddValue(global_index_i, local_rhs[test_loc_i + test_local_offset]);
246
247 for (unsigned int trial_f = 0; trial_f < trial_functions_dofs_data.dofs_data.size(); ++trial_f)
248 {
249 const auto &trial_dofs_data = trial_functions_dofs_data.dofs_data.at(trial_f).get();
250 const auto &trial_global_dofs = trial_dofs_data.CellsGlobalDOFs.at(dimension).at(cell_index);
251 const auto trial_global_offset_DOFs = trial_functions_dofs_data.global_offsets_DOFs[trial_f];
252 const auto trial_local_offset = trial_functions_dofs_data.local_offsets[trial_f];
253 const auto trial_global_offset_Strongs = trial_functions_dofs_data.global_offsets_Strongs[trial_f];
254
255 for (unsigned int trial_loc_j = 0; trial_loc_j < trial_global_dofs.size(); ++trial_loc_j)
256 {
257 const auto &global_dof_j = trial_global_dofs.at(trial_loc_j);
258 const auto &local_dof_j =
259 trial_dofs_data.CellsDOFs.at(global_dof_j.Dimension).at(global_dof_j.CellIndex).at(global_dof_j.DOFIndex);
260
261 const unsigned int global_index_j = local_dof_j.Global_Index;
262 const double loc_A_element = local_lhs(test_loc_i + test_local_offset, trial_loc_j + trial_local_offset);
263
264 switch (local_dof_j.Type)
265 {
267 global_lhs_Strongs.Triplet(global_index_i, global_index_j + trial_global_offset_Strongs, loc_A_element);
268 break;
270 global_lhs_DOFs.Triplet(global_index_i, global_index_j + trial_global_offset_DOFs, loc_A_element);
271 break;
272 default:
273 throw std::runtime_error("Unknown DOF Type");
274 }
275 }
276 }
277 }
278 }
279}
280// ***************************************************************************
281template <unsigned int dimension, typename local_lhs_type, typename global_lhs_type>
282void assemble_local_matrix_to_global_matrix(const unsigned int cell_index,
283 const local_matrix_to_global_matrix_dofs_data &test_functions_dofs_data,
284 const local_matrix_to_global_matrix_dofs_data &trial_functions_dofs_data,
285 const local_lhs_type &local_lhs,
286 global_lhs_type &global_lhs_DOFs,
287 global_lhs_type &global_lhs_Strongs)
288{
289 for (size_t test_f = 0; test_f < test_functions_dofs_data.dofs_data.size(); ++test_f)
290 {
291 const auto &test_dofs_data = test_functions_dofs_data.dofs_data.at(test_f).get();
292 const auto &test_global_dofs = test_dofs_data.CellsGlobalDOFs.at(dimension).at(cell_index);
293 const auto test_global_offset_DOFs = test_functions_dofs_data.global_offsets_DOFs[test_f];
294 const auto test_local_offset = test_functions_dofs_data.local_offsets[test_f];
295
296 for (size_t test_loc_i = 0; test_loc_i < test_global_dofs.size(); ++test_loc_i)
297 {
298 const auto global_dof_i = test_global_dofs.at(test_loc_i);
299 const auto local_dof_i =
300 test_dofs_data.CellsDOFs.at(global_dof_i.Dimension).at(global_dof_i.CellIndex).at(global_dof_i.DOFIndex);
301
302 switch (local_dof_i.Type)
303 {
305 continue;
307 break;
308 default:
309 throw std::runtime_error("Unknown DOF Type");
310 }
311
312 const unsigned int global_index_i = local_dof_i.Global_Index + test_global_offset_DOFs;
313
314 for (unsigned int trial_f = 0; trial_f < trial_functions_dofs_data.dofs_data.size(); ++trial_f)
315 {
316 const auto &trial_dofs_data = trial_functions_dofs_data.dofs_data.at(trial_f).get();
317 const auto &trial_global_dofs = trial_dofs_data.CellsGlobalDOFs.at(dimension).at(cell_index);
318 const auto trial_global_offset_DOFs = trial_functions_dofs_data.global_offsets_DOFs[trial_f];
319 const auto trial_local_offset = trial_functions_dofs_data.local_offsets[trial_f];
320 const auto trial_global_offset_Strongs = trial_functions_dofs_data.global_offsets_Strongs[trial_f];
321
322 for (unsigned int trial_loc_j = 0; trial_loc_j < trial_global_dofs.size(); ++trial_loc_j)
323 {
324 const auto &global_dof_j = trial_global_dofs.at(trial_loc_j);
325 const auto &local_dof_j =
326 trial_dofs_data.CellsDOFs.at(global_dof_j.Dimension).at(global_dof_j.CellIndex).at(global_dof_j.DOFIndex);
327
328 const unsigned int global_index_j = local_dof_j.Global_Index;
329 const double loc_A_element = local_lhs(test_loc_i + test_local_offset, trial_loc_j + trial_local_offset);
330
331 switch (local_dof_j.Type)
332 {
334 global_lhs_Strongs.Triplet(global_index_i, global_index_j + trial_global_offset_Strongs, loc_A_element);
335 break;
337 global_lhs_DOFs.Triplet(global_index_i, global_index_j + trial_global_offset_DOFs, loc_A_element);
338 break;
339 default:
340 throw std::runtime_error("Unknown DOF Type");
341 }
342 }
343 }
344 }
345 }
346}
347// ***************************************************************************
348template <unsigned int dimension, typename local_rhs_type, typename global_rhs_type>
349void assemble_local_matrix_to_global_matrix(const unsigned int cell_index,
350 const local_matrix_to_global_matrix_dofs_data &test_functions_dofs_data,
351 const local_rhs_type &local_rhs,
352 global_rhs_type &global_rhs)
353{
354 for (size_t test_f = 0; test_f < test_functions_dofs_data.dofs_data.size(); ++test_f)
355 {
356 const auto &test_dofs_data = test_functions_dofs_data.dofs_data.at(test_f).get();
357 const auto &test_global_dofs = test_dofs_data.CellsGlobalDOFs.at(dimension).at(cell_index);
358 const auto test_global_offset_DOFs = test_functions_dofs_data.global_offsets_DOFs[test_f];
359 const auto test_local_offset = test_functions_dofs_data.local_offsets[test_f];
360
361 for (size_t test_loc_i = 0; test_loc_i < test_global_dofs.size(); ++test_loc_i)
362 {
363 const auto global_dof_i = test_global_dofs.at(test_loc_i);
364 const auto local_dof_i =
365 test_dofs_data.CellsDOFs.at(global_dof_i.Dimension).at(global_dof_i.CellIndex).at(global_dof_i.DOFIndex);
366
367 switch (local_dof_i.Type)
368 {
370 continue;
372 break;
373 default:
374 throw std::runtime_error("Unknown DOF Type");
375 }
376
377 const unsigned int global_index_i = local_dof_i.Global_Index + test_global_offset_DOFs;
378
379 global_rhs.AddValue(global_index_i, local_rhs[test_loc_i + test_local_offset]);
380 }
381 }
382}
383
384// ***************************************************************************
385} // namespace Assembler_Utilities
386} // namespace PDETools
387} // namespace Polydim
388
389#endif
Polydim::PDETools::Assembler_Utilities::local_count_dofs_data local_count_dofs(const unsigned int cell_index, const std::vector< std::reference_wrapper< const Polydim::PDETools::DOFs::DOFsManager::DOFsData > > &dofs_data)
Definition Assembler_Utilities.hpp:74
void assemble_local_matrix_to_global_matrix(const unsigned int cell_index, const local_matrix_to_global_matrix_dofs_data &test_functions_dofs_data, const local_matrix_to_global_matrix_dofs_data &trial_functions_dofs_data, const local_lhs_type &local_lhs, const local_rhs_type &local_rhs, global_lhs_type &global_lhs_DOFs, global_lhs_type &global_lhs_Strongs, global_rhs_type &global_rhs)
Definition Assembler_Utilities.hpp:211
Eigen::VectorXd global_solution_to_local_solution(const unsigned int cell_index, const std::vector< std::reference_wrapper< const Polydim::PDETools::DOFs::DOFsManager::DOFsData > > &dofs_data, const unsigned int &num_local_DOFs, const std::vector< size_t > &local_offsets_DOFs, const std::vector< size_t > &global_offsets_DOFs, const std::vector< size_t > &global_offsets_Strongs, const global_solution_type &global_solution_DOFs, const global_solution_type &global_solution_Strongs)
Definition Assembler_Utilities.hpp:119
Polydim::PDETools::Assembler_Utilities::count_dofs_data count_dofs(const std::vector< Polydim::PDETools::DOFs::DOFsManager::DOFsData > &dofs_data)
Definition Assembler_Utilities.hpp:47
Definition FEM_MCC_2D_LocalSpace.cpp:17
Definition Assembler_Utilities.hpp:33
unsigned int num_total_dofs
Definition Assembler_Utilities.hpp:35
unsigned int num_total_strong
Definition Assembler_Utilities.hpp:36
std::vector< size_t > offsets_DOFs
Definition Assembler_Utilities.hpp:37
unsigned int num_total_boundary_dofs
Definition Assembler_Utilities.hpp:34
std::vector< size_t > offsets_Strongs
Definition Assembler_Utilities.hpp:38
std::vector< size_t > offsets_DOFs
Definition Assembler_Utilities.hpp:44
unsigned int num_total_dofs
Definition Assembler_Utilities.hpp:43
std::vector< size_t > global_offsets_DOFs
Definition Assembler_Utilities.hpp:28
std::vector< size_t > global_offsets_Strongs
Definition Assembler_Utilities.hpp:29
std::vector< std::reference_wrapper< const Polydim::PDETools::DOFs::DOFsManager::DOFsData > > dofs_data
Definition Assembler_Utilities.hpp:26
std::vector< size_t > local_offsets
Definition Assembler_Utilities.hpp:27
Full DOF numbering produced for a mesh: local/global indices and aggregate counts.
Definition DOFsManager.hpp:116