Hermes2D  3.0
discrete_problem_thread_assembler.cpp
1 // This file is part of Hermes2D.
2 //
3 // Hermes2D is free software: you can redistribute it and/or modify
4 // it under the terms of the GNU General Public License as published by
5 // the Free Software Foundation, either version 2 of the License, or
6 // (at your option) any later version.
7 //
8 // Hermes2D is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 // GNU General Public License for more details.
12 //
13 // You should have received a copy of the GNU General Public License
14 // along with Hermes2D. If not, see <http://www.gnu.org/licenses/>.
15 
16 #include "discrete_problem/discrete_problem_thread_assembler.h"
17 #include "discrete_problem/discrete_problem_selective_assembler.h"
18 #include "shapeset/precalc.h"
19 #include "function/solution.h"
20 #include "weakform/weakform.h"
21 #include "function/exact_solution.h"
22 
23 namespace Hermes
24 {
25  namespace Hermes2D
26  {
27  template<typename Scalar>
28  DiscreteProblemThreadAssembler<Scalar>::DiscreteProblemThreadAssembler(DiscreteProblemSelectiveAssembler<Scalar>* selectiveAssembler, bool nonlinear) :
29  pss(nullptr), refmaps(nullptr), u_ext(nullptr),
30  selectiveAssembler(selectiveAssembler), integrationOrderCalculator(selectiveAssembler),
31  ext_funcs(nullptr), ext_funcs_allocated_size(0), ext_funcs_local(nullptr), ext_funcs_local_allocated_size(0),
32  funcs_wf_initialized(false), funcs_space_initialized(false), spaces_size(0), nonlinear(nonlinear), reusable_DOFs(nullptr), reusable_Dirichlet(nullptr)
33  {
34  // Init the memory pool - if PJLIB is linked, it will do the magic, if not, it will initialize the pointer to null.
35  this->init_funcs_memory_pool();
36  }
37 
38  template<typename Scalar>
39  DiscreteProblemThreadAssembler<Scalar>::~DiscreteProblemThreadAssembler()
40  {
41  this->free();
42 #ifdef WITH_PJLIB
43  pj_pool_release(this->FuncMemoryPool);
44 #endif
45  }
46 
47  template<typename Scalar>
48  void DiscreteProblemThreadAssembler<Scalar>::init_spaces(const std::vector<SpaceSharedPtr<Scalar> > spaces)
49  {
50  this->free_spaces();
51 
52  bool reinit_funcs = this->spaces_size != spaces.size();
53  this->spaces_size = spaces.size();
54 
55  if (reinit_funcs)
56  {
57  this->deinit_funcs_space();
58  this->init_funcs_space();
59  }
60 
61  pss = malloc_with_check<PrecalcShapesetAssembling*>(spaces_size);
62  refmaps = malloc_with_check<RefMap*>(spaces_size);
63 
64  for (unsigned int j = 0; j < spaces_size; j++)
65  {
66  pss[j] = new PrecalcShapesetAssembling(spaces[j]->shapeset);
67  refmaps[j] = new RefMap();
68  refmaps[j]->set_quad_2d(&g_quad_2d_std);
69  }
70  }
71 
72  template<typename Scalar>
73  void DiscreteProblemThreadAssembler<Scalar>::set_weak_formulation(WeakFormSharedPtr<Scalar> wf_)
74  {
75  this->deinit_funcs_wf();
76  this->free_weak_formulation();
77  this->wf = WeakFormSharedPtr<Scalar>(wf_->clone());
78  this->wf->cloneMembers(wf_);
79  this->init_funcs_wf();
80  }
81 
82  template<typename Scalar>
83  void DiscreteProblemThreadAssembler<Scalar>::init_u_ext(const std::vector<SpaceSharedPtr<Scalar> > spaces, Solution<Scalar>** u_ext_sln)
84  {
85  assert(this->spaces_size == spaces.size() && this->pss);
86 
87  free_u_ext();
88  u_ext = malloc_with_check<Solution<Scalar>*>(spaces_size);
89 
90  for (unsigned int j = 0; j < spaces_size; j++)
91  {
92  if (u_ext_sln)
93  {
94  u_ext[j] = new Solution<Scalar>(spaces[j]->get_mesh());
95  u_ext[j]->copy(u_ext_sln[j]);
96  }
97  else
98  {
99  if (spaces[j]->get_shapeset()->get_num_components() == 1)
100  u_ext[j] = new ZeroSolution<Scalar>(spaces[j]->get_mesh());
101  else
102  u_ext[j] = new ZeroSolutionVector<Scalar>(spaces[j]->get_mesh());
103  }
104  }
105 
106  this->integrationOrderCalculator.u_ext = this->u_ext;
107  }
108 
109  template<typename Scalar>
110  void DiscreteProblemThreadAssembler<Scalar>::init_assembling(Solution<Scalar>** u_ext_sln, const std::vector<SpaceSharedPtr<Scalar> >& spaces, bool add_dirichlet_lift_)
111  {
112  // Basic settings.
113  this->add_dirichlet_lift = add_dirichlet_lift_;
114 
115  // Transformables setup.
116  fns.clear();
117  // - precalc shapesets.
118  for (unsigned j = 0; j < this->spaces_size; j++)
119  {
120  fns.push_back(pss[j]);
121  pss[j]->set_quad_2d(&g_quad_2d_std);
122  }
123  // - wf->ext.
124  for (unsigned j = 0; j < this->wf->ext.size(); j++)
125  {
126  fns.push_back(this->wf->ext[j].get());
127  this->wf->ext[j]->set_quad_2d(&g_quad_2d_std);
128  }
129  // - forms->ext.
130  for (unsigned int form_i = 0; form_i < this->wf->get_forms().size(); form_i++)
131  {
132  Form<Scalar>* form = this->wf->get_forms()[form_i];
133  for (unsigned int ext_i = 0; ext_i < form->ext.size(); ext_i++)
134  {
135  if (form->ext[ext_i])
136  {
137  fns.push_back(form->ext[ext_i].get());
138  form->ext[ext_i]->set_quad_2d(&g_quad_2d_std);
139  }
140  }
141  }
142  // - u_ext.
143  if (this->nonlinear)
144  {
145  init_u_ext(spaces, u_ext_sln);
146  for (unsigned j = 0; j < this->wf->get_neq(); j++)
147  {
148  fns.push_back(u_ext[j]);
149  u_ext[j]->set_quad_2d(&g_quad_2d_std);
150  }
151  }
152 
153  // Process markers.
154  this->wf->processFormMarkers(spaces);
155  }
156 
157  template<typename Scalar>
158  void DiscreteProblemThreadAssembler<Scalar>::init_funcs_memory_pool()
159  {
160 #ifdef WITH_PJLIB
161  pj_thread_desc rtpdesc;
162  pj_thread_t *thread;
163  pj_bzero(rtpdesc, sizeof(rtpdesc));
164  if (!pj_thread_is_registered())
165  pj_thread_register(NULL, rtpdesc, &thread);
166 
167  this->FuncMemoryPool = pj_pool_create(&Hermes2DMemoryPoolCache.factory, // the factory
168  "pool-DiscreteProblemThreadAssembler", // pool's name
169  sizeof(Func<Scalar>) * H2D_MAX_LOCAL_BASIS_SIZE * (H2D_MAX_NUMBER_EDGES + 1), // initial size
170  sizeof(Func<Scalar>) * H2D_MAX_LOCAL_BASIS_SIZE * H2D_MAX_NUMBER_EDGES, // increment size
171  NULL);
172 #else
173  this->FuncMemoryPool = nullptr;
174 #endif
175  }
176 
177  template<typename Scalar>
178  void DiscreteProblemThreadAssembler<Scalar>::init_funcs_space()
179  {
180  this->funcs_space_initialized = true;
181 
182  // Basis & test fns, u_ext funcs.
183  for (unsigned short space_i = 0; space_i < this->spaces_size; space_i++)
184  {
185  for (unsigned int j = 0; j < H2D_MAX_LOCAL_BASIS_SIZE; j++)
186  this->funcs[space_i][j] = preallocate_fn<double>(this->FuncMemoryPool);
187 
188  for (int edge_i = 0; edge_i < H2D_MAX_NUMBER_EDGES; edge_i++)
189  for (unsigned int j = 0; j < H2D_MAX_LOCAL_BASIS_SIZE; j++)
190  this->funcsSurface[edge_i][space_i][j] = preallocate_fn<double>(this->FuncMemoryPool);
191 
192  if (this->nonlinear)
193  this->u_ext_funcs[space_i] = preallocate_fn<Scalar>(this->FuncMemoryPool);
194  }
195  }
196 
197  template<typename Scalar>
198  void DiscreteProblemThreadAssembler<Scalar>::init_funcs_wf()
199  {
200  this->funcs_wf_initialized = true;
201 
202  // Reallocation of wf-(nonlocal-) ext funcs.
203  int ext_size = this->wf->ext.size();
204  int u_ext_fns_size = this->wf->u_ext_fn.size();
205 
206  if (ext_size + u_ext_fns_size > ext_funcs_allocated_size)
207  {
208  ext_funcs_allocated_size = ext_size + u_ext_fns_size;
209  ext_funcs = realloc_with_check<DiscreteProblemThreadAssembler<Scalar>, Func<Scalar>*>(ext_funcs, ext_funcs_allocated_size, this);
210  }
211 
212  // Initializaton of wf-(nonlocal-)ext funcs
213  if (ext_size > 0 || u_ext_fns_size > 0)
214  {
215  for (int ext_i = 0; ext_i < u_ext_fns_size; ext_i++)
216  this->ext_funcs[ext_i] = preallocate_fn<Scalar>(this->FuncMemoryPool);
217 
218  for (int ext_i = 0; ext_i < ext_size; ext_i++)
219  this->ext_funcs[u_ext_fns_size + ext_i] = preallocate_fn<Scalar>(this->FuncMemoryPool);
220  }
221 
222  // Calculating local sizes.
223  unsigned short local_ext_size = 0;
224  unsigned short local_u_ext_fns_size = 0;
225  for (unsigned short form_i = 0; form_i < this->wf->forms.size(); form_i++)
226  {
227  if (this->wf->forms[form_i]->ext.size() > local_ext_size)
228  local_ext_size = this->wf->forms[form_i]->ext.size();
229 
230  if (this->wf->forms[form_i]->u_ext_fn.size() > local_u_ext_fns_size)
231  local_u_ext_fns_size = this->wf->forms[form_i]->u_ext_fn.size();
232  }
233 
234  if (local_ext_size > 0 || local_u_ext_fns_size > 0)
235  {
236  // Reallocation of form-(local-) ext funcs.
237  if (local_ext_size + local_u_ext_fns_size > ext_funcs_local_allocated_size)
238  {
239  ext_funcs_local_allocated_size = local_ext_size + local_u_ext_fns_size;
240  ext_funcs_local = realloc_with_check<DiscreteProblemThreadAssembler<Scalar>, Func<Scalar>*>(ext_funcs_local, ext_funcs_local_allocated_size, this);
241  }
242 
243  // Initializaton of form-(local-)ext funcs
244  for (int ext_i = 0; ext_i < local_u_ext_fns_size; ext_i++)
245  this->ext_funcs_local[ext_i] = preallocate_fn<Scalar>(this->FuncMemoryPool);
246 
247  for (int ext_i = 0; ext_i < local_ext_size; ext_i++)
248  this->ext_funcs_local[local_u_ext_fns_size + ext_i] = preallocate_fn<Scalar>(this->FuncMemoryPool);
249  }
250  }
251 
252  template<typename Scalar>
253  void DiscreteProblemThreadAssembler<Scalar>::deinit_funcs_space()
254  {
255  if (!funcs_space_initialized)
256  return;
257  else
258  funcs_space_initialized = false;
259 
260  for (unsigned short space_i = 0; space_i < this->spaces_size; space_i++)
261  {
262  // Test functions
263  for (unsigned int j = 0; j < H2D_MAX_LOCAL_BASIS_SIZE; j++)
264  delete this->funcs[space_i][j];
265 
266  // Test functions - surface
267  for (int edge_i = 0; edge_i < H2D_MAX_NUMBER_EDGES; edge_i++)
268  {
269  for (unsigned int j = 0; j < H2D_MAX_LOCAL_BASIS_SIZE; j++)
270  delete this->funcsSurface[edge_i][space_i][j];
271  }
272 
273  // UExt
274  if (this->nonlinear)
275  delete this->u_ext_funcs[space_i];
276  }
277  }
278 
279  template<typename Scalar>
280  void DiscreteProblemThreadAssembler<Scalar>::deinit_funcs_wf()
281  {
282  if (!funcs_wf_initialized)
283  return;
284  else
285  funcs_wf_initialized = false;
286 
287  // Ext
288  int ext_size = this->wf->ext.size();
289  int u_ext_fns_size = this->wf->u_ext_fn.size();
290  if (ext_size > 0 || u_ext_fns_size > 0)
291  {
292  for (int ext_i = 0; ext_i < u_ext_fns_size; ext_i++)
293  delete this->ext_funcs[ext_i];
294 
295  for (int ext_i = 0; ext_i < ext_size; ext_i++)
296  delete this->ext_funcs[u_ext_fns_size + ext_i];
297  }
298 
299  // Ext - local
300  unsigned short local_ext_size = 0;
301  unsigned short local_u_ext_fns_size = 0;
302  for (unsigned short form_i = 0; form_i < this->wf->forms.size(); form_i++)
303  {
304  if (this->wf->forms[form_i]->ext.size() > local_ext_size)
305  local_ext_size = this->wf->forms[form_i]->ext.size();
306 
307  if (this->wf->forms[form_i]->u_ext_fn.size() > local_u_ext_fns_size)
308  local_u_ext_fns_size = this->wf->forms[form_i]->u_ext_fn.size();
309  }
310 
311  if (local_ext_size > 0 || local_u_ext_fns_size > 0)
312  {
313  for (int ext_i = 0; ext_i < local_u_ext_fns_size; ext_i++)
314  delete this->ext_funcs_local[ext_i];
315 
316  for (int ext_i = 0; ext_i < local_ext_size; ext_i++)
317  delete this->ext_funcs_local[local_u_ext_fns_size + ext_i];
318  }
319  }
320 
321  template<typename Scalar>
322  void DiscreteProblemThreadAssembler<Scalar>::deinit_funcs()
323  {
324  this->deinit_funcs_wf();
325  this->deinit_funcs_space();
326  }
327 
328  template<typename Scalar>
329  void DiscreteProblemThreadAssembler<Scalar>::init_assembling_one_state(const std::vector<SpaceSharedPtr<Scalar> >& spaces, Traverse::State* current_state_)
330  {
331  current_state = current_state_;
332  this->integrationOrderCalculator.current_state = this->current_state;
333 
334  // Active elements.
335  for (unsigned short j = 0; j < fns.size(); j++)
336  {
337  if (current_state->e[j])
338  {
339  fns[j]->set_active_element(current_state->e[j]);
340  fns[j]->set_transform(current_state->sub_idx[j]);
341  }
342  }
343 
344  // Assembly lists & refmaps.
345  for (int j = 0; j < this->spaces_size; j++)
346  {
347  if (current_state->e[j])
348  {
349  spaces[j]->get_element_assembly_list(current_state->e[j], &als[j]);
350  refmaps[j]->set_active_element(current_state->e[j]);
351  refmaps[j]->force_transform(pss[j]->get_transform(), pss[j]->get_ctm());
352  rep_refmap = refmaps[j];
353  }
354  }
355 
356  // Boundary assembly lists
357  if (current_state->isBnd && !(this->wf->mfsurf.empty() && this->wf->vfsurf.empty()))
358  {
359  for (int j = 0; j < this->spaces_size; j++)
360  {
361  if (current_state->e[j])
362  {
363  for (int k = 0; k < current_state->rep->nvert; k++)
364  {
365  if (current_state->bnd[k])
366  spaces[j]->get_boundary_assembly_list(current_state->e[j], k, &alsSurface[k][j]);
367  }
368  }
369  }
370  }
371 
372  // Volumetric integration order.
373  this->order = this->integrationOrderCalculator.calculate_order(spaces, this->refmaps, this->wf);
374 
375  // Init the variables (funcs, geometry, ...)
376  this->init_calculation_variables();
377  }
378 
379  template<typename Scalar>
380  void DiscreteProblemThreadAssembler<Scalar>::init_calculation_variables()
381  {
382  for (unsigned short space_i = 0; space_i < this->spaces_size; space_i++)
383  {
384  if (current_state->e[space_i] == nullptr)
385  continue;
386 
387  for (unsigned int j = 0; j < this->als[space_i].cnt; j++)
388  {
389  pss[space_i]->set_active_shape(this->als[space_i].idx[j]);
390  init_fn_preallocated(this->funcs[space_i][j], pss[space_i], refmaps[space_i], this->order);
391  }
392  }
393 
394  this->n_quadrature_points = init_geometry_points_allocated(this->rep_refmap, this->order, this->geometry, this->jacobian_x_weights);
395 
396  if (current_state->isBnd && (this->wf->mfsurf.size() > 0 || this->wf->vfsurf.size() > 0))
397  {
398  int order_local = this->order;
399 
400  for (unsigned int edge_i = 0; edge_i < this->current_state->rep->nvert; edge_i++)
401  {
402  if (!current_state->bnd[edge_i])
403  continue;
404 
405  this->n_quadrature_pointsSurface[edge_i] = init_surface_geometry_points_allocated(this->rep_refmap, this->order, edge_i, current_state->rep->marker, this->geometrySurface[edge_i], this->jacobian_x_weightsSurface[edge_i]);
406  this->orderSurface[edge_i] = this->order;
407  this->order = order_local;
408 
409  for (unsigned short space_i = 0; space_i < this->spaces_size; space_i++)
410  {
411  if (!current_state->e[space_i])
412  continue;
413 
414  for (unsigned int j = 0; j < this->alsSurface[edge_i][space_i].cnt; j++)
415  {
416  pss[space_i]->set_active_shape(this->alsSurface[edge_i][space_i].idx[j]);
417  init_fn_preallocated(this->funcsSurface[edge_i][space_i][j], pss[space_i], refmaps[space_i], this->orderSurface[edge_i]);
418  }
419  }
420  }
421  }
422  }
423 
424  template<typename Scalar>
425  void DiscreteProblemThreadAssembler<Scalar>::deinit_calculation_variables()
426  {
427  }
428 
429  template<typename Scalar>
430  void DiscreteProblemThreadAssembler<Scalar>::init_u_ext_values(int order)
431  {
432  if (this->nonlinear)
433  {
434  for (int i = 0; i < spaces_size; i++)
435  {
436  if (u_ext[i]->get_active_element())
437  init_fn_preallocated(u_ext_funcs[i], u_ext[i], order);
438  }
439  }
440  }
441 
442  template<typename Scalar>
443  template<typename Geom>
444  void DiscreteProblemThreadAssembler<Scalar>::init_ext_values(Func<Scalar>** target_array, std::vector<MeshFunctionSharedPtr<Scalar> >& ext, std::vector<UExtFunctionSharedPtr<Scalar> >& u_ext_fns, int order, Func<Scalar>** u_ext_func, Geom* geometry)
445  {
446  int ext_size = ext.size();
447  int u_ext_fns_size = u_ext_fns.size();
448 
449  if (ext_size > 0 || u_ext_fns_size > 0)
450  {
451  for (int ext_i = 0; ext_i < ext_size; ext_i++)
452  {
453  if (ext[ext_i])
454  {
455  if (ext[ext_i]->get_active_element())
456  init_fn_preallocated(target_array[u_ext_fns_size + ext_i], ext[ext_i].get(), order);
457  }
458  }
459 
460  for (int ext_i = 0; ext_i < u_ext_fns_size; ext_i++)
461  {
462  if (u_ext_fns[ext_i])
463  init_fn_preallocated(target_array[ext_i], u_ext_fns[ext_i].get(), target_array, u_ext_func, order, geometry, current_state->rep->get_mode());
464  }
465  }
466 
467  if (this->rungeKutta)
468  {
469  for (unsigned short ext_i = 0; ext_i < ext.size(); ext_i++)
470  u_ext_func[ext_i]->add(target_array[ext.size() - this->RK_original_spaces_count + ext_i]);
471  }
472  }
473 
474  template<typename Scalar>
475  void DiscreteProblemThreadAssembler<Scalar>::assemble_one_state()
476  {
477  // init - u_ext_func
478  this->init_u_ext_values(this->order);
479 
480  // init - ext
481  this->init_ext_values(this->ext_funcs, this->wf->ext, this->wf->u_ext_fn, this->order, this->u_ext_funcs, &this->geometry);
482 
483  if (this->current_mat || this->add_dirichlet_lift)
484  {
485  for (unsigned short current_mfvol_i = 0; current_mfvol_i < this->wf->mfvol.size(); current_mfvol_i++)
486  {
487  if (!selectiveAssembler->form_to_be_assembled(this->wf->mfvol[current_mfvol_i], current_state))
488  continue;
489 
490  int form_i = this->wf->mfvol[current_mfvol_i]->i;
491  int form_j = this->wf->mfvol[current_mfvol_i]->j;
492 
493  this->assemble_matrix_form(this->wf->mfvol[current_mfvol_i], order, funcs[form_j], funcs[form_i], &als[form_i], &als[form_j], n_quadrature_points, &geometry, jacobian_x_weights);
494  }
495  }
496  if (this->current_rhs)
497  {
498  for (unsigned short current_vfvol_i = 0; current_vfvol_i < this->wf->vfvol.size(); current_vfvol_i++)
499  {
500  if (!selectiveAssembler->form_to_be_assembled(this->wf->vfvol[current_vfvol_i], current_state))
501  continue;
502 
503  int form_i = this->wf->vfvol[current_vfvol_i]->i;
504 
505  this->assemble_vector_form(this->wf->vfvol[current_vfvol_i], order, funcs[form_i], &als[form_i], n_quadrature_points, &geometry, jacobian_x_weights);
506  }
507  }
508 
509  // Assemble surface integrals now: loop through surfaces of the element.
510  if (current_state->isBnd && (this->wf->mfsurf.size() > 0 || this->wf->vfsurf.size() > 0))
511  {
512  for (unsigned char isurf = 0; isurf < current_state->rep->nvert; isurf++)
513  {
514  if (!current_state->bnd[isurf])
515  continue;
516 
517  current_state->isurf = isurf;
518 
519  // Edge-wise parameters for WeakForm.
520  this->wf->set_active_edge_state(current_state->e, isurf);
521 
522  // init - u_ext_func
523  this->init_u_ext_values(this->orderSurface[isurf]);
524 
525  // init - ext
526  this->init_ext_values(this->ext_funcs, this->wf->ext, this->wf->u_ext_fn, this->orderSurface[isurf], this->u_ext_funcs, &this->geometrySurface[isurf]);
527 
528  if (this->current_mat || this->add_dirichlet_lift)
529  {
530  for (unsigned short current_mfsurf_i = 0; current_mfsurf_i < this->wf->mfsurf.size(); current_mfsurf_i++)
531  {
532  if (!selectiveAssembler->form_to_be_assembled(this->wf->mfsurf[current_mfsurf_i], current_state))
533  continue;
534 
535  int form_i = this->wf->mfsurf[current_mfsurf_i]->i;
536  int form_j = this->wf->mfsurf[current_mfsurf_i]->j;
537 
538  this->assemble_matrix_form(this->wf->mfsurf[current_mfsurf_i], orderSurface[isurf], funcsSurface[isurf][form_j], funcsSurface[isurf][form_i],
539  &alsSurface[isurf][form_i], &alsSurface[isurf][form_j], n_quadrature_pointsSurface[isurf], &geometrySurface[isurf], jacobian_x_weightsSurface[isurf]);
540  }
541  }
542 
543  if (this->current_rhs)
544  {
545  for (unsigned short current_vfsurf_i = 0; current_vfsurf_i < this->wf->vfsurf.size(); current_vfsurf_i++)
546  {
547  if (!selectiveAssembler->form_to_be_assembled(this->wf->vfsurf[current_vfsurf_i], current_state))
548  continue;
549 
550  int form_i = this->wf->vfsurf[current_vfsurf_i]->i;
551 
552  this->assemble_vector_form(this->wf->vfsurf[current_vfsurf_i], orderSurface[isurf], funcsSurface[isurf][form_i], &alsSurface[isurf][form_i],
553  n_quadrature_pointsSurface[isurf], &geometrySurface[isurf], jacobian_x_weightsSurface[isurf]);
554  }
555  }
556  }
557  }
558  }
559 
560  template<typename Scalar>
561  template<typename MatrixFormType, typename Geom>
562  void DiscreteProblemThreadAssembler<Scalar>::assemble_matrix_form(MatrixFormType* form, int order, Func<double>** base_fns, Func<double>** test_fns,
563  AsmList<Scalar>* current_als_i, AsmList<Scalar>* current_als_j, int n_quadrature_points, Geom* geometry, double* jacobian_x_weights)
564  {
565  bool surface_form = (dynamic_cast<MatrixFormVol<Scalar>*>(form) == nullptr);
566 
567  double block_scaling_coefficient = this->block_scaling_coeff(form);
568 
569  bool tra = (form->i != form->j) && (form->sym != 0);
570  bool sym = (form->i == form->j) && (form->sym == 1);
571 
572  Func<Scalar>** ext_local = this->ext_funcs;
573  // If the user supplied custom ext functions for this form.
574  if (form->ext.size() > 0 || form->u_ext_fn.size() > 0)
575  {
576  this->init_ext_values(this->ext_funcs_local, form->ext, (form->u_ext_fn.size() > 0 ? form->u_ext_fn : this->wf->u_ext_fn), order, this->u_ext_funcs, geometry);
577  ext_local = this->ext_funcs_local;
578  }
579 
580  // Account for the previous time level solution previously inserted at the back of ext.
581  Func<Scalar>** u_ext_local = this->u_ext_funcs;
582  if (this->rungeKutta)
583  u_ext_local += form->u_ext_offset;
584 
585  // Actual form-specific calculation.
586  for (unsigned int i = 0; i < current_als_i->cnt; i++)
587  {
588  if (current_als_i->dof[i] < 0 || std::abs(current_als_i->coef[i]) < Hermes::HermesSqrtEpsilon)
589  continue;
590 
591  if ((!tra || surface_form) && current_als_i->dof[i] < 0)
592  continue;
593 
594  for (unsigned int j = 0; j < current_als_j->cnt; j++)
595  {
596  if (current_als_j->dof[j] >= 0 && this->reusable_DOFs && *this->reusable_DOFs)
597  {
598  if ((*this->reusable_DOFs)[current_als_j->dof[j]] && (*this->reusable_DOFs)[current_als_i->dof[i]])
599  {
600  unsigned short local_matrix_index_array = i * H2D_MAX_LOCAL_BASIS_SIZE + j;
601  local_stiffness_matrix[local_matrix_index_array] = 0.;
602  if (sym || tra)
603  {
604  unsigned short local_matrix_index_array_transposed = j * H2D_MAX_LOCAL_BASIS_SIZE + i;
605  local_stiffness_matrix[local_matrix_index_array_transposed] = 0.;
606  }
607  continue;
608  }
609  }
610 
611  if (current_als_j->dof[j] < 0 && this->reusable_Dirichlet && *this->reusable_Dirichlet && (*this->reusable_Dirichlet)[form->j])
612  {
613  continue;
614  }
615 
616  // Skip symmetric values that do not contribute to Dirichlet lift.
617  if (sym && j < i && current_als_j->dof[j] >= 0)
618  continue;
619 
620  // Skip anything that does not contribute to Dirichlet in the case of just rhs assembling.
621  if (current_als_j->dof[j] >= 0 && !this->current_mat)
622  continue;
623 
624  if (std::abs(current_als_j->coef[j]) < Hermes::HermesEpsilon)
625  continue;
626 
627  Func<double>* u = base_fns[j];
628  Func<double>* v = test_fns[i];
629 
630  Scalar val = block_scaling_coefficient * form->value(n_quadrature_points, jacobian_x_weights, u_ext_local, u, v, geometry, ext_local) * form->scaling_factor * current_als_j->coef[j] * current_als_i->coef[i];
631 
632  if (current_als_j->dof[j] >= 0)
633  {
634  unsigned short local_matrix_index_array = i * H2D_MAX_LOCAL_BASIS_SIZE + j;
635 
636  if (surface_form)
637  local_stiffness_matrix[local_matrix_index_array] = 0.5 * val;
638  else
639  local_stiffness_matrix[local_matrix_index_array] = val;
640 
641  if (sym)
642  {
643  unsigned short local_matrix_index_array_transposed = j * H2D_MAX_LOCAL_BASIS_SIZE + i;
644  local_stiffness_matrix[local_matrix_index_array_transposed] = local_stiffness_matrix[local_matrix_index_array];
645  }
646  }
647  else if (this->add_dirichlet_lift && this->current_rhs)
648  {
649  this->dirichlet_lift_rhs->add(current_als_i->dof[i], -val);
650  }
651  }
652  }
653 
654  // Insert the local stiffness matrix into the global one.
655  if (this->current_mat)
656  this->current_mat->add(current_als_i->cnt, current_als_j->cnt, local_stiffness_matrix, current_als_i->dof, current_als_j->dof, H2D_MAX_LOCAL_BASIS_SIZE);
657 
658  // Insert also the off-diagonal (anti-)symmetric block, if required.
659  if (tra)
660  {
661  if (form->sym < 0)
662  change_sign(local_stiffness_matrix, current_als_i->cnt, current_als_j->cnt, H2D_MAX_LOCAL_BASIS_SIZE);
663  transpose(local_stiffness_matrix, current_als_i->cnt, current_als_j->cnt, H2D_MAX_LOCAL_BASIS_SIZE);
664 
665  if (this->current_mat)
666  this->current_mat->add(current_als_j->cnt, current_als_i->cnt, local_stiffness_matrix, current_als_j->dof, current_als_i->dof, H2D_MAX_LOCAL_BASIS_SIZE);
667 
668  if (this->add_dirichlet_lift && this->current_rhs)
669  {
670  for (unsigned int j = 0; j < current_als_i->cnt; j++)
671  {
672  if (current_als_i->dof[j] < 0)
673  {
674  for (unsigned int i = 0; i < current_als_j->cnt; i++)
675  {
676  if (current_als_j->dof[i] >= 0)
677  {
678  int local_matrix_index_array = i * H2D_MAX_LOCAL_BASIS_SIZE + j;
679  this->dirichlet_lift_rhs->add(current_als_j->dof[i], -local_stiffness_matrix[local_matrix_index_array]);
680  }
681  }
682  }
683  }
684  }
685  }
686  }
687 
688  template<typename Scalar>
689  template<typename VectorFormType, typename Geom>
690  void DiscreteProblemThreadAssembler<Scalar>::assemble_vector_form(VectorFormType* form, int order, Func<double>** test_fns,
691  AsmList<Scalar>* current_als_i, int n_quadrature_points, Geom* geometry, double* jacobian_x_weights)
692  {
693  bool surface_form = (dynamic_cast<VectorFormVol<Scalar>*>(form) == nullptr);
694 
695  Func<Scalar>** ext_local = this->ext_funcs;
696  // If the user supplied custom ext functions for this form.
697  if (form->ext.size() > 0 || form->u_ext_fn.size() > 0)
698  {
699  this->init_ext_values(this->ext_funcs_local, form->ext, (form->u_ext_fn.size() > 0 ? form->u_ext_fn : this->wf->u_ext_fn), order, this->u_ext_funcs, geometry);
700  ext_local = this->ext_funcs_local;
701  }
702 
703  // Account for the previous time level solution previously inserted at the back of ext.
704  Func<Scalar>** u_ext_local = this->u_ext_funcs;
705  if (this->rungeKutta)
706  u_ext_local += form->u_ext_offset;
707 
708  // Actual form-specific calculation.
709  for (unsigned int i = 0; i < current_als_i->cnt; i++)
710  {
711  if (current_als_i->dof[i] < 0)
712  continue;
713 
714  if (this->reusable_DOFs && *this->reusable_DOFs)
715  {
716  if ((*this->reusable_DOFs)[current_als_i->dof[i]])
717  continue;
718  }
719 
720  // Is this necessary, i.e. is there a coefficient smaller than Hermes::HermesSqrtEpsilon?
721  if (std::abs(current_als_i->coef[i]) < Hermes::HermesSqrtEpsilon)
722  continue;
723 
724  Func<double>* v = test_fns[i];
725 
726  Scalar val;
727  if (surface_form)
728  val = 0.5 * form->value(n_quadrature_points, jacobian_x_weights, u_ext_local, v, geometry, ext_local) * form->scaling_factor * current_als_i->coef[i];
729  else
730  val = form->value(n_quadrature_points, jacobian_x_weights, u_ext_local, v, geometry, ext_local) * form->scaling_factor * current_als_i->coef[i];
731 
732  this->current_rhs->add(current_als_i->dof[i], val);
733  }
734  }
735 
736  template<typename Scalar>
737  void DiscreteProblemThreadAssembler<Scalar>::deinit_assembling_one_state()
738  {
739  this->deinit_calculation_variables();
740  }
741 
742  template<typename Scalar>
743  void DiscreteProblemThreadAssembler<Scalar>::deinit_assembling()
744  {
745  }
746 
747  template<typename Scalar>
749  {
750  this->deinit_funcs();
751  this->free_spaces();
752  this->free_weak_formulation();
753  this->free_u_ext();
754 
755  free_with_check(ext_funcs, true);
756  free_with_check(ext_funcs_local, true);
757  }
758 
759  template<typename Scalar>
761  {
762  if (!this->pss)
763  return;
764 
765  for (unsigned int j = 0; j < spaces_size; j++)
766  delete pss[j];
767  free_with_check(pss);
768 
769  for (unsigned int j = 0; j < spaces_size; j++)
770  delete refmaps[j];
771  free_with_check(refmaps);
772  }
773 
774  template<typename Scalar>
775  void DiscreteProblemThreadAssembler<Scalar>::free_weak_formulation()
776  {
777  if (this->wf)
778  this->wf->free_ext();
779  }
780 
781  template<typename Scalar>
782  void DiscreteProblemThreadAssembler<Scalar>::free_u_ext()
783  {
784  if (u_ext)
785  {
786  for (unsigned int j = 0; j < spaces_size; j++)
787  delete u_ext[j];
788  free_with_check(u_ext);
789  }
790  }
791 
792  template class HERMES_API DiscreteProblemThreadAssembler < double > ;
793  template class HERMES_API DiscreteProblemThreadAssembler < std::complex<double> > ;
794  }
795 }
Definition: adapt.h:24
#define H2D_MAX_NUMBER_EDGES
A maximum number of edges of an element.
Definition: global.h:31
HERMES_API void init_fn_preallocated(Func< double > *u, PrecalcShapeset *fu, RefMap *rm, const int order)
Init the shape function for the evaluation of the volumetric/surface integral (transformation of valu...
Definition: forms.cpp:497
This class is a one-thread (non-DG) assembly worker.
HERMES_API unsigned char init_geometry_points_allocated(RefMap **reference_mapping, unsigned short reference_mapping_count, int order, GeomVol< double > &geometry, double *jacobian_x_weights)