Hermes2D  3.0
weakform.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 "global.h"
17 #include "api2d.h"
18 #include "weakform.h"
19 #include "matrix.h"
20 #include "forms.h"
21 #include "shapeset_l2_all.h"
22 #include "shapeset_hc_all.h"
23 #include "shapeset_hd_all.h"
24 #include "shapeset_h1_all.h"
25 #include "algebra/dense_matrix_operations.h"
26 
28 
29 namespace Hermes
30 {
31  namespace Hermes2D
32  {
33  template<typename Scalar>
34  WeakFormSharedPtr<Scalar>::WeakFormSharedPtr(Hermes::Hermes2D::WeakForm<Scalar>* ptr) : std::tr1::shared_ptr<Hermes::Hermes2D::WeakForm<Scalar> >(ptr)
35  {
36  }
37 
38  template<typename Scalar>
39  WeakFormSharedPtr<Scalar>::WeakFormSharedPtr(const WeakFormSharedPtr& other) : std::tr1::shared_ptr<Hermes::Hermes2D::WeakForm<Scalar> >(other)
40  {
41  }
42 
43  template<typename Scalar>
44  void WeakFormSharedPtr<Scalar>::operator=(const WeakFormSharedPtr& other)
45  {
46  std::tr1::shared_ptr<Hermes::Hermes2D::WeakForm<Scalar> >::operator=(other);
47  }
48 
49  static bool warned_nonOverride = false;
50 
51  // This is to be used by weak forms specifying numerical flux through interior edges.
52  // Forms with this identifier will receive DiscontinuousFunc representations of shape
53  // and ext. functions, which they may query for values on either side of given interface.
54  static const std::string H2D_DG_INNER_EDGE = "-1234567";
55 
56  template<typename Scalar>
57  WeakForm<Scalar>::WeakForm(unsigned int neq, bool mat_free) : Hermes::Mixins::Loggable(false)
58  {
59  this->neq = neq;
60  this->original_neq = neq;
61  this->is_matfree = mat_free;
62  }
63 
64  template<typename Scalar>
66  {
67  this->ext.clear();
68  }
69 
70  template<typename Scalar>
72  {
73  for (unsigned int i = 0; i < this->forms.size(); i++)
74  delete get_forms()[i];
75  delete_all();
76  }
77 
78  template<typename Scalar>
80  {
81  if (!warned_nonOverride)
82 #pragma omp critical (warning_weakform_nonOverride)
83  {
84  if (!warned_nonOverride)
85  this->warn("Using default WeakForm<Scalar>::clone, if you have any dynamically created data in your WeakForm constructor, you need to overload this method!");
86  warned_nonOverride = true;
87  }
88  return new WeakForm(*this);
89  }
90 
91  template<typename Scalar>
93  {
94  // Nothing here, supposed to be overriden if necessary.
95  }
96 
97  template<typename Scalar>
98  void WeakForm<Scalar>::set_active_edge_state(Element** e, unsigned char isurf)
99  {
100  // Nothing here, supposed to be overriden if necessary.
101  }
102 
103  template<typename Scalar>
104  void WeakForm<Scalar>::set_active_DG_state(Element** e, unsigned char isurf)
105  {
106  // Nothing here, supposed to be overriden if necessary.
107  }
108 
109  template<typename Scalar>
111  {
112  this->mfvol.clear();
113  this->mfsurf.clear();
114  this->mfDG.clear();
115  this->vfvol.clear();
116  this->vfsurf.clear();
117  this->vfDG.clear();
118  this->forms.clear();
119  this->ext.clear();
120 
121  for (unsigned int i = 0; i < other_wf->forms.size(); i++)
122  {
123  if (dynamic_cast<MatrixFormVol<Scalar>*>(other_wf->forms[i]) != nullptr)
124  {
125  MatrixFormVol<Scalar>* form = (dynamic_cast<MatrixFormVol<Scalar>*>(other_wf->forms[i]))->clone();
126  form->copy_base(other_wf->forms[i]);
127  this->forms.push_back(form);
128  this->mfvol.push_back(dynamic_cast<MatrixFormVol<Scalar>*>(this->forms.back()));
129  continue;
130  }
131  if (dynamic_cast<VectorFormVol<Scalar>*>(other_wf->forms[i]) != nullptr)
132  {
133  VectorFormVol<Scalar>* form = (dynamic_cast<VectorFormVol<Scalar>*>(other_wf->forms[i]))->clone();
134  form->copy_base(other_wf->forms[i]);
135  this->forms.push_back(form);
136  this->vfvol.push_back(dynamic_cast<VectorFormVol<Scalar>*>(this->forms.back()));
137  continue;
138  }
139  if (dynamic_cast<MatrixFormSurf<Scalar>*>(other_wf->forms[i]) != nullptr)
140  {
141  MatrixFormSurf<Scalar>* form = (dynamic_cast<MatrixFormSurf<Scalar>*>(other_wf->forms[i]))->clone();
142  form->copy_base(other_wf->forms[i]);
143  this->forms.push_back(form);
144  this->mfsurf.push_back(dynamic_cast<MatrixFormSurf<Scalar>*>(this->forms.back()));
145  continue;
146  }
147  if (dynamic_cast<VectorFormSurf<Scalar>*>(other_wf->forms[i]) != nullptr)
148  {
149  VectorFormSurf<Scalar>* form = (dynamic_cast<VectorFormSurf<Scalar>*>(other_wf->forms[i]))->clone();
150  form->copy_base(other_wf->forms[i]);
151  this->forms.push_back(form);
152  this->vfsurf.push_back(dynamic_cast<VectorFormSurf<Scalar>*>(this->forms.back()));
153  continue;
154  }
155  if (dynamic_cast<MatrixFormDG<Scalar>*>(other_wf->forms[i]) != nullptr)
156  {
157  MatrixFormDG<Scalar>* form = (dynamic_cast<MatrixFormDG<Scalar>*>(other_wf->forms[i]))->clone();
158  form->copy_base(other_wf->forms[i]);
159  this->forms.push_back(form);
160  this->mfDG.push_back(dynamic_cast<MatrixFormDG<Scalar>*>(this->forms.back()));
161  continue;
162  }
163  if (dynamic_cast<VectorFormDG<Scalar>*>(other_wf->forms[i]) != nullptr)
164  {
165  VectorFormDG<Scalar>* form = (dynamic_cast<VectorFormDG<Scalar>*>(other_wf->forms[i]))->clone();
166  form->copy_base(other_wf->forms[i]);
167  this->forms.push_back(form);
168  this->forms.push_back((dynamic_cast<VectorFormDG<Scalar>*>(other_wf->forms[i]))->clone());
169  this->vfDG.push_back(dynamic_cast<VectorFormDG<Scalar>*>(this->forms.back()));
170  continue;
171  }
172  }
173  for (unsigned int i = 0; i < other_wf->forms.size(); i++)
174  {
175  this->cloneMemberExtFunctions(other_wf->forms[i]->ext, this->forms[i]->ext);
176  other_wf->forms[i]->u_ext_fn = this->forms[i]->u_ext_fn;
177  this->forms[i]->wf = this;
178  }
179  this->cloneMemberExtFunctions(other_wf->ext, this->ext);
180  this->u_ext_fn = other_wf->u_ext_fn;
181  }
182 
183  template<typename Scalar>
184  void WeakForm<Scalar>::cloneMemberExtFunctions(std::vector<MeshFunctionSharedPtr<Scalar> > source_ext, std::vector<MeshFunctionSharedPtr<Scalar> >& cloned_ext)
185  {
186  cloned_ext.clear();
187  for (unsigned int i = 0; i < source_ext.size(); i++)
188  {
189  Solution<Scalar>* originalSln = dynamic_cast<Solution<Scalar>*>(source_ext[i].get());
190  if (originalSln)
191  {
192  Solution<Scalar>* newSln = nullptr;
193  if (originalSln->get_type() == HERMES_SLN)
194  {
195  newSln = new Solution < Scalar > ;
196  newSln->copy(source_ext[i].get());
197  }
198  else
199  newSln = static_cast<Solution<Scalar>*>(originalSln->clone());
200 
201  newSln->set_type(originalSln->get_type());
202  cloned_ext.push_back(newSln);
203  }
204  else
205  cloned_ext.push_back(source_ext[i]->clone());
206  }
207  }
208 
209  template<typename Scalar>
211  {
212  mfvol.clear();
213  mfsurf.clear();
214  mfDG.clear();
215  vfvol.clear();
216  vfsurf.clear();
217  vfDG.clear();
218  forms.clear();
219  };
220 
221  template<typename Scalar>
223  {
224  this->ext.clear();
225  this->ext.push_back(ext);
226  }
227 
228  template<typename Scalar>
230  {
231  this->ext = ext;
232  }
233 
234  template<typename Scalar>
236  {
237  this->u_ext_fn.clear();
238  this->u_ext_fn.push_back(ext);
239  }
240 
241  template<typename Scalar>
243  {
244  this->u_ext_fn = ext;
245  }
246 
247  template<typename Scalar>
248  template<typename FormType>
249  void WeakForm<Scalar>::processFormMarkers(const std::vector<SpaceSharedPtr<Scalar> > spaces, bool surface, std::vector<FormType> forms_to_process)
250  {
251  for (unsigned short form_i = 0; form_i < forms_to_process.size(); form_i++)
252  {
253  Form<Scalar>* form = forms_to_process[form_i];
254  form->areas_internal.clear();
255  for (unsigned short marker_i = 0; marker_i < form->areas.size(); marker_i++)
256  {
257  if (form->areas[marker_i] == HERMES_ANY)
258  {
259  form->assembleEverywhere = true;
260  form->areas_internal.clear();
261  break;
262  }
263 
264  Mesh::MarkersConversion::IntValid marker;
265  if (surface)
266  marker = spaces[form->i]->get_mesh()->get_boundary_markers_conversion().get_internal_marker(form->areas[marker_i]);
267  else
268  marker = spaces[form->i]->get_mesh()->get_element_markers_conversion().get_internal_marker(form->areas[marker_i]);
269 
270  if (marker.valid)
271  form->areas_internal.push_back(marker.marker);
272  else
273  throw Exceptions::Exception("Marker not valid in assembling: %s.", form->areas[marker_i].c_str());
274  }
275  }
276  }
277 
278  template<typename Scalar>
279  void WeakForm<Scalar>::processFormMarkers(const std::vector<SpaceSharedPtr<Scalar> > spaces)
280  {
281  processFormMarkers(spaces, false, this->mfvol);
282  processFormMarkers(spaces, false, this->vfvol);
283  processFormMarkers(spaces, true, this->mfsurf);
284  processFormMarkers(spaces, true, this->vfsurf);
285  }
286 
287  template<typename Scalar>
288  bool WeakForm<Scalar>::is_DG() const
289  {
290  if (this->mfDG.empty() && this->vfDG.empty())
291  return false;
292  return true;
293  }
294 
295  template<typename Scalar>
296  std::vector<MeshFunctionSharedPtr<Scalar> > WeakForm<Scalar>::get_ext() const
297  {
298  return this->ext;
299  }
300 
301  template<typename Scalar>
302  Form<Scalar>::Form(int i) : scaling_factor(1.0), wf(nullptr), assembleEverywhere(false), i(i)
303  {
304  areas.push_back(HERMES_ANY);
305  stage_time = 0.0;
306  }
307 
308  template<typename Scalar>
310  {
311  }
312 
313  template<typename Scalar>
314  void Form<Scalar>::set_current_stage_time(double time)
315  {
316  stage_time = time;
317  }
318 
319  template<typename Scalar>
320  double Form<Scalar>::get_current_stage_time() const
321  {
322  return stage_time;
323  }
324 
325  template<typename Scalar>
327  {
328  areas.clear();
329  areas.push_back(area);
330  }
331  template<typename Scalar>
332  void Form<Scalar>::set_areas(std::vector<std::string> areas)
333  {
334  this->areas = areas;
335  }
336 
337  template<typename Scalar>
338  std::vector<std::string> Form<Scalar>::getAreas() const
339  {
340  return this->areas;
341  }
342 
343  template<typename Scalar>
344  void Form<Scalar>::setScalingFactor(double scalingFactor)
345  {
346  this->scaling_factor = scalingFactor;
347  }
348 
349  template<typename Scalar>
351  {
352  this->ext.clear();
353  this->ext.push_back(ext);
354  }
355 
356  template<typename Scalar>
358  {
359  this->ext = ext;
360  }
361 
362  template<typename Scalar>
364  {
365  this->u_ext_fn.clear();
366  this->u_ext_fn.push_back(ext);
367  }
368 
369  template<typename Scalar>
370  void Form<Scalar>::set_u_ext_fn(std::vector<UExtFunctionSharedPtr<Scalar> > ext)
371  {
372  this->u_ext_fn = ext;
373  }
374 
375  template<typename Scalar>
376  std::vector<MeshFunctionSharedPtr<Scalar> > Form<Scalar>::get_ext() const
377  {
378  return this->ext;
379  }
380 
381  template<typename Scalar>
382  void Form<Scalar>::copy_base(Form<Scalar>* other_form)
383  {
384  this->stage_time = other_form->stage_time;
385  this->scaling_factor = other_form->scaling_factor;
386  this->u_ext_offset = other_form->u_ext_offset;
387  this->previous_iteration_space_index = other_form->previous_iteration_space_index;
388  }
389 
390  template<typename Scalar>
391  MatrixForm<Scalar>::MatrixForm(unsigned int i, unsigned int j) :
392  Form<Scalar>(i), sym(HERMES_NONSYM), j(j)
393  {
395  }
396 
397  template<typename Scalar>
399  {
400  }
401 
402  template<typename Scalar>
403  MatrixFormVol<Scalar>::MatrixFormVol(unsigned int i, unsigned int j) :
404  MatrixForm<Scalar>(i, j)
405  {
406  }
407 
408  template<typename Scalar>
410  {
411  }
412 
413  template<typename Scalar>
414  void MatrixFormVol<Scalar>::setSymFlag(SymFlag sym)
415  {
416  this->sym = sym;
417  }
418 
419  template<typename Scalar>
420  SymFlag MatrixFormVol<Scalar>::getSymFlag() const
421  {
422  return this->sym;
423  }
424 
425  template<typename Scalar>
426  Scalar MatrixFormVol<Scalar>::value(int n, double *wt, Func<Scalar> **u_ext, Func<double> *u, Func<double> *v,
427  GeomVol<double> *e, Func<Scalar> **ext) const
428  {
429  throw Hermes::Exceptions::MethodNotOverridenException("MatrixFormVol<Scalar>::value");
430  return 0.0;
431  }
432 
433  template<typename Scalar>
434  Hermes::Ord MatrixFormVol<Scalar>::ord(int n, double *wt, Func<Hermes::Ord> **u_ext, Func<Hermes::Ord> *u, Func<Hermes::Ord> *v,
435  GeomVol<Hermes::Ord> *e, Func<Ord> **ext) const
436  {
437  throw Hermes::Exceptions::MethodNotOverridenException("MatrixFormVol<Scalar>::ord");
438  return Hermes::Ord();
439  }
440 
441  template<typename Scalar>
442  MatrixFormVol<Scalar>* MatrixFormVol<Scalar>::clone() const
443  {
444  throw Hermes::Exceptions::MethodNotOverridenException("MatrixFormVol<Scalar>::clone()");
445  return nullptr;
446  }
447 
448  template<typename Scalar>
449  MatrixFormSurf<Scalar>::MatrixFormSurf(unsigned int i, unsigned int j) :
450  MatrixForm<Scalar>(i, j)
451  {
452  }
453 
454  template<typename Scalar>
455  Scalar MatrixFormSurf<Scalar>::value(int n, double *wt, Func<Scalar> **u_ext, Func<double> *u, Func<double> *v,
456  GeomSurf<double> *e, Func<Scalar> **ext) const
457  {
458  throw Hermes::Exceptions::MethodNotOverridenException("MatrixFormSurf<Scalar>::value");
459  return 0.0;
460  }
461 
462  template<typename Scalar>
463  Hermes::Ord MatrixFormSurf<Scalar>::ord(int n, double *wt, Func<Hermes::Ord> **u_ext, Func<Hermes::Ord> *u, Func<Hermes::Ord> *v,
464  GeomSurf<Hermes::Ord> *e, Func<Ord> **ext) const
465  {
466  throw Hermes::Exceptions::MethodNotOverridenException("MatrixFormSurf<Scalar>::ord");
467  return Hermes::Ord();
468  }
469 
470  template<typename Scalar>
471  MatrixFormSurf<Scalar>* MatrixFormSurf<Scalar>::clone() const
472  {
473  throw Hermes::Exceptions::MethodNotOverridenException("MatrixFormSurf<Scalar>::clone()");
474  return nullptr;
475  }
476 
477  template<typename Scalar>
478  MatrixFormSurf<Scalar>::~MatrixFormSurf()
479  {
480  }
481 
482  template<typename Scalar>
483  MatrixFormDG<Scalar>::MatrixFormDG(unsigned int i, unsigned int j) :
484  Form<Scalar>(i), j(j)
485  {
487  this->set_area(H2D_DG_INNER_EDGE);
488  }
489 
490  template<typename Scalar>
492  {
493  }
494 
495  template<typename Scalar>
496  Scalar MatrixFormDG<Scalar>::value(int n, double *wt, DiscontinuousFunc<Scalar> **u_ext, DiscontinuousFunc<double> *u, DiscontinuousFunc<double> *v,
497  InterfaceGeom<double> *e, DiscontinuousFunc<Scalar> **ext) const
498  {
499  throw Hermes::Exceptions::MethodNotOverridenException("MatrixFormDG<Scalar>::value");
500  return 0.0;
501  }
502 
503  template<typename Scalar>
504  Hermes::Ord MatrixFormDG<Scalar>::ord(int n, double *wt, DiscontinuousFunc<Hermes::Ord> **u_ext, DiscontinuousFunc<Hermes::Ord> *u, DiscontinuousFunc<Hermes::Ord> *v,
505  InterfaceGeom<Hermes::Ord> *e, DiscontinuousFunc<Ord> **ext) const
506  {
507  throw Hermes::Exceptions::MethodNotOverridenException("MatrixFormDG<Scalar>::ord");
508  return Hermes::Ord();
509  }
510 
511  template<typename Scalar>
512  MatrixFormDG<Scalar>* MatrixFormDG<Scalar>::clone() const
513  {
514  throw Hermes::Exceptions::MethodNotOverridenException("MatrixFormDG<Scalar>::clone()");
515  return nullptr;
516  }
517 
518  template<typename Scalar>
520  Form<Scalar>(i)
521  {
523  }
524 
525  template<typename Scalar>
527  {
528  }
529 
530  template<typename Scalar>
532  VectorForm<Scalar>(i)
533  {
534  }
535 
536  template<typename Scalar>
538  {
539  }
540 
541  template<typename Scalar>
542  Scalar VectorFormVol<Scalar>::value(int n, double *wt, Func<Scalar> **u_ext, Func<double> *v,
543  GeomVol<double> *e, Func<Scalar> **ext) const
544  {
545  throw Hermes::Exceptions::MethodNotOverridenException("VectorFormVol<Scalar>::value");
546  return 0.0;
547  }
548 
549  template<typename Scalar>
550  Hermes::Ord VectorFormVol<Scalar>::ord(int n, double *wt, Func<Hermes::Ord> **u_ext, Func<Hermes::Ord> *v,
551  GeomVol<Hermes::Ord> *e, Func<Ord> **ext) const
552  {
553  throw Hermes::Exceptions::MethodNotOverridenException("VectorFormVol<Scalar>::ord");
554  return Hermes::Ord();
555  }
556 
557  template<typename Scalar>
558  VectorFormVol<Scalar>* VectorFormVol<Scalar>::clone() const
559  {
560  throw Hermes::Exceptions::MethodNotOverridenException("VectorFormVol<Scalar>::clone()");
561  return nullptr;
562  }
563 
564  template<typename Scalar>
566  VectorForm<Scalar>(i)
567  {
568  }
569 
570  template<typename Scalar>
571  Scalar VectorFormSurf<Scalar>::value(int n, double *wt, Func<Scalar> **u_ext, Func<double> *v,
572  GeomSurf<double> *e, Func<Scalar> **ext) const
573  {
574  throw Hermes::Exceptions::MethodNotOverridenException("VectorFormSurf<Scalar>::value");
575  return 0.0;
576  }
577 
578  template<typename Scalar>
579  Hermes::Ord VectorFormSurf<Scalar>::ord(int n, double *wt, Func<Hermes::Ord> **u_ext, Func<Hermes::Ord> *v,
580  GeomSurf<Hermes::Ord> *e, Func<Ord> **ext) const
581  {
582  throw Hermes::Exceptions::MethodNotOverridenException("VectorFormSurf<Scalar>::ord");
583  return Hermes::Ord();
584  }
585 
586  template<typename Scalar>
587  VectorFormSurf<Scalar>::~VectorFormSurf()
588  {
589  }
590 
591  template<typename Scalar>
592  VectorFormSurf<Scalar>* VectorFormSurf<Scalar>::clone() const
593  {
594  throw Hermes::Exceptions::MethodNotOverridenException("VectorFormSurf<Scalar>::clone()");
595  return nullptr;
596  }
597 
598  template<typename Scalar>
600  Form<Scalar>(i)
601  {
602  this->set_area(H2D_DG_INNER_EDGE);
603  }
604 
605  template<typename Scalar>
607  {
608  }
609 
610  template<typename Scalar>
611  Scalar VectorFormDG<Scalar>::value(int n, double *wt, DiscontinuousFunc<Scalar> **u_ext, Func<double> *v,
612  InterfaceGeom<double> *e, DiscontinuousFunc<Scalar> **ext) const
613  {
614  throw Hermes::Exceptions::MethodNotOverridenException("VectorFormDG<Scalar>::value");
615  return 0.0;
616  }
617 
618  template<typename Scalar>
619  Hermes::Ord VectorFormDG<Scalar>::ord(int n, double *wt, DiscontinuousFunc<Hermes::Ord> **u_ext, Func<Hermes::Ord> *v,
620  InterfaceGeom<Hermes::Ord> *e, DiscontinuousFunc<Ord> **ext) const
621  {
622  throw Hermes::Exceptions::MethodNotOverridenException("VectorFormDG<Scalar>::ord");
623  return Hermes::Ord();
624  }
625 
626  template<typename Scalar>
627  VectorFormDG<Scalar>* VectorFormDG<Scalar>::clone() const
628  {
629  throw Hermes::Exceptions::MethodNotOverridenException("VectorFormDG<Scalar>::clone()");
630  return nullptr;
631  }
632 
633  template<typename Scalar>
635  {
636  this->wf = wf;
637  if (wf->original_neq != wf->neq)
638  this->previous_iteration_space_index %= wf->original_neq;
639  }
640 
641  template<typename Scalar>
643  {
644  if (form->i >= neq || form->j >= neq)
645  throw Hermes::Exceptions::Exception("Invalid equation number.");
646  if (form->sym < -1 || form->sym > 1)
647  throw Hermes::Exceptions::Exception("\"sym\" must be -1, 0 or 1.");
648  if (form->sym < 0 && form->i == form->j)
649  throw Hermes::Exceptions::Exception("Only off-diagonal forms can be antisymmetric.");
650  if (mfvol.size() > 100)
651  {
652  this->warn("Large number of forms (> 100). Is this the intent?");
653  }
654 
655  form->set_weakform(this);
656  mfvol.push_back(form);
657  forms.push_back(form);
658  }
659 
660  template<typename Scalar>
662  {
663  if (form->i >= neq || form->j >= neq)
664  throw Hermes::Exceptions::Exception("Invalid equation number.");
665 
666  form->set_weakform(this);
667  mfsurf.push_back(form);
668  forms.push_back(form);
669  }
670 
671  template<typename Scalar>
673  {
674  if (form->i >= neq || form->j >= neq)
675  throw Hermes::Exceptions::Exception("Invalid equation number.");
676 
677  form->set_weakform(this);
678  mfDG.push_back(form);
679  forms.push_back(form);
680  }
681 
682  template<typename Scalar>
684  {
685  if (form->i >= neq)
686  throw Hermes::Exceptions::Exception("Invalid equation number.");
687  form->set_weakform(this);
688  vfvol.push_back(form);
689  forms.push_back(form);
690  }
691 
692  template<typename Scalar>
694  {
695  if (form->i >= neq)
696  throw Hermes::Exceptions::Exception("Invalid equation number.");
697 
698  form->set_weakform(this);
699  vfsurf.push_back(form);
700  forms.push_back(form);
701  }
702 
703  template<typename Scalar>
705  {
706  if (form->i >= neq)
707  throw Hermes::Exceptions::Exception("Invalid equation number.");
708 
709  form->set_weakform(this);
710  vfDG.push_back(form);
711  forms.push_back(form);
712  }
713 
714  template<typename Scalar>
715  std::vector<Form<Scalar> *> WeakForm<Scalar>::get_forms() const
716  {
717  return forms;
718  }
719 
720  template<typename Scalar>
721  std::vector<MatrixFormVol<Scalar> *> WeakForm<Scalar>::get_mfvol() const
722  {
723  return mfvol;
724  }
725  template<typename Scalar>
726  std::vector<MatrixFormSurf<Scalar> *> WeakForm<Scalar>::get_mfsurf() const
727  {
728  return mfsurf;
729  }
730  template<typename Scalar>
731  std::vector<MatrixFormDG<Scalar> *> WeakForm<Scalar>::get_mfDG() const
732  {
733  return mfDG;
734  }
735  template<typename Scalar>
736  std::vector<VectorFormVol<Scalar> *> WeakForm<Scalar>::get_vfvol() const
737  {
738  return vfvol;
739  }
740  template<typename Scalar>
741  std::vector<VectorFormSurf<Scalar> *> WeakForm<Scalar>::get_vfsurf() const
742  {
743  return vfsurf;
744  }
745  template<typename Scalar>
746  std::vector<VectorFormDG<Scalar> *> WeakForm<Scalar>::get_vfDG() const
747  {
748  return vfDG;
749  }
750 
751  template<typename Scalar>
752  bool** WeakForm<Scalar>::get_blocks(bool force_diagonal_blocks) const
753  {
754  bool** blocks = new_matrix<bool>(neq, neq);
755  for (unsigned int i = 0; i < neq; i++)
756  {
757  for (unsigned int j = 0; j < neq; j++)
758  blocks[i][j] = false;
759  if (force_diagonal_blocks)
760  blocks[i][i] = true;
761  }
762  for (unsigned i = 0; i < mfvol.size(); i++)
763  {
764  if (fabs(mfvol[i]->scaling_factor) > Hermes::HermesSqrtEpsilon)
765  blocks[mfvol[i]->i][mfvol[i]->j] = true;
766  if (mfvol[i]->sym)
767  if (fabs(mfvol[i]->scaling_factor) > Hermes::HermesSqrtEpsilon)
768  blocks[mfvol[i]->j][mfvol[i]->i] = true;
769  }
770  for (unsigned i = 0; i < mfsurf.size(); i++)
771  {
772  if (fabs(mfsurf[i]->scaling_factor) > Hermes::HermesSqrtEpsilon)
773  blocks[mfsurf[i]->i][mfsurf[i]->j] = true;
774  }
775  for (unsigned i = 0; i < mfDG.size(); i++)
776  {
777  if (fabs(mfDG[i]->scaling_factor) > Hermes::HermesSqrtEpsilon)
778  blocks[mfDG[i]->i][mfDG[i]->j] = true;
779  }
780 
781  return blocks;
782  }
783 
784  template<typename Scalar>
786  {
787  current_time = time;
788  }
789 
790  template<typename Scalar>
792  {
793  return current_time;
794  }
795 
796  template<typename Scalar>
798  {
799  current_time_step = time_step;
800  }
801 
802  template<typename Scalar>
804  {
805  return current_time_step;
806  }
807 
808  template class HERMES_API WeakForm < double > ;
809  template class HERMES_API WeakForm < std::complex<double> > ;
810  template class HERMES_API Form < double > ;
811  template class HERMES_API Form < std::complex<double> > ;
812  template class HERMES_API MatrixForm < double > ;
813  template class HERMES_API MatrixForm < std::complex<double> > ;
814  template class HERMES_API MatrixFormVol < double > ;
815  template class HERMES_API MatrixFormVol < std::complex<double> > ;
816  template class HERMES_API MatrixFormSurf < double > ;
817  template class HERMES_API MatrixFormSurf < std::complex<double> > ;
818  template class HERMES_API MatrixFormDG < double > ;
819  template class HERMES_API MatrixFormDG < std::complex<double> > ;
820  template class HERMES_API VectorForm < double > ;
821  template class HERMES_API VectorForm < std::complex<double> > ;
822  template class HERMES_API VectorFormVol < double > ;
823  template class HERMES_API VectorFormVol < std::complex<double> > ;
824  template class HERMES_API VectorFormSurf < double > ;
825  template class HERMES_API VectorFormSurf < std::complex<double> > ;
826  template class HERMES_API VectorFormDG < double > ;
827  template class HERMES_API VectorFormDG < std::complex<double> > ;
828 
829  template class HERMES_API WeakFormSharedPtr < double > ;
830  template class HERMES_API WeakFormSharedPtr < std::complex<double> > ;
831  }
832 }
VectorFormSurf(unsigned int i)
Constructor with coordinates.
Definition: weakform.cpp:565
VectorFormDG(unsigned int i)
Constructor with coordinates.
Definition: weakform.cpp:599
void add_matrix_form_surf(MatrixFormSurf< Scalar > *mfs)
Adds surface matrix form.
Definition: weakform.cpp:661
MatrixFormVol(unsigned int i, unsigned int j)
Constructor with coordinates.
Definition: weakform.cpp:403
void setScalingFactor(double scalingFactor)
scaling factor
Definition: weakform.cpp:344
unsigned int original_neq
Original number of equations in case this is a Runge-Kutta enlarged system.
Definition: weakform.h:213
Definition: adapt.h:24
Abstract, base class for vector form - i.e. a single integral in the linear form on the right hand si...
Definition: weakform.h:439
void set_area(std::string area)
Definition: weakform.cpp:326
Stores one element of a mesh.
Definition: element.h:107
void set_ext(MeshFunctionSharedPtr< Scalar > ext)
Definition: weakform.cpp:222
void set_ext(MeshFunctionSharedPtr< Scalar > ext)
Definition: weakform.cpp:350
void add_vector_form_DG(VectorFormDG< Scalar > *vfDG)
Adds DG vector form.
Definition: weakform.cpp:704
virtual void set_active_state(Element **e)
Definition: weakform.cpp:92
unsigned int neq
Number of equations.
Definition: weakform.h:210
MatrixFormDG(unsigned int i, unsigned int j)
Constructor with coordinates.
Definition: weakform.cpp:483
std::vector< int > areas_internal
Internal - this structure is being filled anew with every assembling.
Definition: weakform.h:310
Used to pass the instances of Space around.
Definition: space.h:34
void add_matrix_form_DG(MatrixFormDG< Scalar > *mfDG)
Adds DG matrix form.
Definition: weakform.cpp:672
virtual WeakForm * clone() const
Cloning.
Definition: weakform.cpp:79
Common definitions for Hermes2D.
::xsd::cxx::tree::time< char, simple_type > time
C++ type corresponding to the time XML Schema built-in type.
virtual void set_active_DG_state(Element **e, unsigned char isurf)
Definition: weakform.cpp:104
void add_vector_form(VectorFormVol< Scalar > *vfv)
Adds volumetric vector form.
Definition: weakform.cpp:683
VectorFormVol(unsigned int i)
Constructor with coordinates.
Definition: weakform.cpp:531
Abstract, base class for vector Surface form - i.e. VectorForm, where the integration is with respect...
Definition: weakform.h:48
MatrixForm(unsigned int i, unsigned int j)
Constructor with coordinates.
Definition: weakform.cpp:391
Abstract, base class for vector DG form - i.e. linear Form, where the integration is with respect to ...
Definition: weakform.h:50
MatrixFormSurf(unsigned int i, unsigned int j)
Constructor with coordinates.
Definition: weakform.cpp:449
Abstract, base class for any form - i.e. integral in the weak formulation of (a system of) PDE By de...
Definition: weakform.h:43
Abstract, base class for matrix Surface form - i.e. MatrixForm, where the integration is with respect...
Definition: weakform.h:47
SymFlag
Bilinear form symmetry flag, see WeakForm::add_matrix_form.
Definition: global.h:156
Represents the weak formulation of a PDE problem.
Definition: global.h:86
void add_vector_form_surf(VectorFormSurf< Scalar > *vfs)
Adds surface vector form.
Definition: weakform.cpp:693
std::vector< Form< Scalar > * > get_forms() const
Internal.
Definition: weakform.cpp:715
virtual double get_current_time() const
Definition: weakform.cpp:791
Abstract, base class for matrix form - i.e. a single integral in the bilinear form on the left hand s...
Definition: weakform.h:357
virtual ~WeakForm()
Destructor.
Definition: weakform.cpp:71
Used to pass the instances of WeakForm around.
Definition: weakform.h:55
virtual double get_current_time_step() const
Definition: weakform.cpp:803
Calculated function values (from the class Function) on an element for assembling.
Definition: forms.h:214
void set_current_time_step(double time_step)
Definition: weakform.cpp:797
::xsd::cxx::tree::string< char, simple_type > string
C++ type corresponding to the string XML Schema built-in type.
void add_matrix_form(MatrixFormVol< Scalar > *mfv)
Adds volumetric matrix form.
Definition: weakform.cpp:642
Abstract, base class for matrix Volumetric form - i.e. MatrixForm, where the integration is with resp...
Definition: weakform.h:45
WeakForm(unsigned int neq=1, bool mat_free=false)
Definition: weakform.cpp:57
unsigned int previous_iteration_space_index
Definition: weakform.h:325
void set_weakform(WeakForm< Scalar > *wf)
Set pointer to a WeakForm + handling of internal data.
Definition: weakform.cpp:634
std::vector< std::string > areas
Markers of the areas where this form will be assembled.
Definition: weakform.h:307
Abstract, base class for matrix DG form - i.e. bilinear form, where the integration is with respect t...
Definition: weakform.h:49
virtual void set_active_edge_state(Element **e, unsigned char isurf)
Definition: weakform.cpp:98
std::vector< MeshFunctionSharedPtr< Scalar > > get_ext() const
Definition: weakform.cpp:296
void set_current_time(double time)
Definition: weakform.cpp:785
void delete_all()
Deletes all volumetric and surface forms.
Definition: weakform.cpp:210
VectorForm(unsigned int i)
Constructor with coordinates.
Definition: weakform.cpp:519
Form(int i=0)
Constructor with coordinates.
Definition: weakform.cpp:302
void set_u_ext_fn(UExtFunctionSharedPtr< Scalar > ext)
Definition: weakform.cpp:235
Abstract, base class for vector Volumetric form - i.e. VectorForm, where the integration is with resp...
Definition: weakform.h:46