Hermes2D  2.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 using namespace Hermes::Algebra::DenseMatrixOperations;
26 namespace Hermes
27 {
28  namespace Hermes2D
29  {
30 
34  static const std::string H2D_DG_INNER_EDGE = "-1234567";
35 
36  template<typename Scalar>
37  WeakForm<Scalar>::WeakForm(unsigned int neq, bool mat_free) : Hermes::Mixins::Loggable(true), warned_nonOverride(false)
38  {
39  this->neq = neq;
40  this->is_matfree = mat_free;
41  }
42 
43  template<typename Scalar>
44  void WeakForm<Scalar>::free_ext()
45  {
46  for(unsigned int i = 0; i < this->ext.size(); i++)
47  delete this->ext[i];
48  for(unsigned int i = 0; i < this->forms.size(); i++)
49  for(unsigned int j = 0; j < get_forms()[i]->ext.size(); j++)
50  delete get_forms()[i]->ext[j];
51  }
52 
53  template<typename Scalar>
54  WeakForm<Scalar>::~WeakForm()
55  {
56  for(unsigned int i = 0; i < this->forms.size(); i++)
57  delete get_forms()[i];
58  delete_all();
59  }
60 
61  template<typename Scalar>
62  WeakForm<Scalar>* WeakForm<Scalar>::clone() const
63  {
64  if(!this->warned_nonOverride)
65 #pragma omp critical (warning_weakform_nonOverride)
66  {
67  if(!this->warned_nonOverride)
68  this->warn("Using default WeakForm<Scalar>::clone, if you have any dynamically created data in your WeakForm constructor, you need to overload this method!");
69  const_cast<WeakForm<Scalar>*>(this)->warned_nonOverride = true;
70  }
71  return new WeakForm(*this);
72  }
73 
74  template<typename Scalar>
75  void WeakForm<Scalar>::cloneMembers(const WeakForm<Scalar>* otherWf)
76  {
77  this->mfvol.clear();
78  this->mfsurf.clear();
79  this->mfDG.clear();
80  this->vfvol.clear();
81  this->vfsurf.clear();
82  this->vfDG.clear();
83  this->forms.clear();
84  this->ext.clear();
85 
86  for(unsigned int i = 0; i < otherWf->forms.size(); i++)
87  {
88  if(dynamic_cast<MatrixFormVol<Scalar>*>(otherWf->forms[i]) != NULL)
89  this->forms.push_back((dynamic_cast<MatrixFormVol<Scalar>*>(otherWf->forms[i]))->clone());
90  if(dynamic_cast<MatrixFormSurf<Scalar>*>(otherWf->forms[i]) != NULL)
91  this->forms.push_back((dynamic_cast<MatrixFormSurf<Scalar>*>(otherWf->forms[i]))->clone());
92  if(dynamic_cast<MatrixFormDG<Scalar>*>(otherWf->forms[i]) != NULL)
93  this->forms.push_back((dynamic_cast<MatrixFormDG<Scalar>*>(otherWf->forms[i]))->clone());
94 
95  if(dynamic_cast<VectorFormVol<Scalar>*>(otherWf->forms[i]) != NULL)
96  this->forms.push_back((dynamic_cast<VectorFormVol<Scalar>*>(otherWf->forms[i]))->clone());
97  if(dynamic_cast<VectorFormSurf<Scalar>*>(otherWf->forms[i]) != NULL)
98  this->forms.push_back((dynamic_cast<VectorFormSurf<Scalar>*>(otherWf->forms[i]))->clone());
99  if(dynamic_cast<VectorFormDG<Scalar>*>(otherWf->forms[i]) != NULL)
100  this->forms.push_back((dynamic_cast<VectorFormDG<Scalar>*>(otherWf->forms[i]))->clone());
101 
102  Hermes::vector<MeshFunction<Scalar>*> newExt;
103  for(unsigned int ext_i = 0; ext_i < otherWf->forms[i]->ext.size(); ext_i++)
104  newExt.push_back(otherWf->forms[i]->ext[ext_i]->clone());
105  this->forms.back()->set_ext(newExt);
106  this->forms.back()->wf = this;
107 
108  if(dynamic_cast<MatrixFormVol<Scalar>*>(otherWf->forms[i]) != NULL)
109  this->mfvol.push_back(dynamic_cast<MatrixFormVol<Scalar>*>(this->forms.back()));
110  if(dynamic_cast<MatrixFormSurf<Scalar>*>(otherWf->forms[i]) != NULL)
111  this->mfsurf.push_back(dynamic_cast<MatrixFormSurf<Scalar>*>(this->forms.back()));
112  if(dynamic_cast<MatrixFormDG<Scalar>*>(otherWf->forms[i]) != NULL)
113  this->mfDG.push_back(dynamic_cast<MatrixFormDG<Scalar>*>(this->forms.back()));
114 
115  if(dynamic_cast<VectorFormVol<Scalar>*>(otherWf->forms[i]) != NULL)
116  this->vfvol.push_back(dynamic_cast<VectorFormVol<Scalar>*>(this->forms.back()));
117  if(dynamic_cast<VectorFormSurf<Scalar>*>(otherWf->forms[i]) != NULL)
118  this->vfsurf.push_back(dynamic_cast<VectorFormSurf<Scalar>*>(this->forms.back()));
119  if(dynamic_cast<VectorFormDG<Scalar>*>(otherWf->forms[i]) != NULL)
120  this->vfDG.push_back(dynamic_cast<VectorFormDG<Scalar>*>(this->forms.back()));
121  }
122  for(unsigned int i = 0; i < otherWf->ext.size(); i++)
123  this->ext.push_back(otherWf->ext[i]->clone());
124  }
125 
126  template<typename Scalar>
127  void WeakForm<Scalar>::delete_all()
128  {
129  mfvol.clear();
130  mfsurf.clear();
131  mfDG.clear();
132  vfvol.clear();
133  vfsurf.clear();
134  vfDG.clear();
135  forms.clear();
136  };
137 
138  template<typename Scalar>
139  void WeakForm<Scalar>::set_ext(MeshFunction<Scalar>* ext)
140  {
141  this->ext.clear();
142  this->ext.push_back(ext);
143  }
144 
145  template<typename Scalar>
146  void WeakForm<Scalar>::set_ext(Hermes::vector<MeshFunction<Scalar>*> ext)
147  {
148  this->ext = ext;
149  }
150 
151  template<typename Scalar>
152  Hermes::vector<MeshFunction<Scalar>*> WeakForm<Scalar>::get_ext() const
153  {
154  return this->ext;
155  }
156 
157  template<typename Scalar>
158  Form<Scalar>::Form() : scaling_factor(1.0), u_ext_offset(0), wf(NULL)
159  {
160  areas.push_back(HERMES_ANY);
161  stage_time = 0.0;
162  }
163 
164  template<typename Scalar>
165  Form<Scalar>::~Form()
166  {
167  }
168 
169  template<typename Scalar>
170  void Form<Scalar>::set_current_stage_time(double time)
171  {
172  stage_time = time;
173  }
174 
175  template<typename Scalar>
176  double Form<Scalar>::get_current_stage_time() const
177  {
178  return stage_time;
179  }
180 
181  template<typename Scalar>
182  void Form<Scalar>::set_area(std::string area)
183  {
184  areas.clear();
185  areas.push_back(area);
186  }
187  template<typename Scalar>
188  void Form<Scalar>::set_areas(Hermes::vector<std::string> areas)
189  {
190  this->areas = areas;
191  }
192 
193  template<typename Scalar>
194  Hermes::vector<std::string> Form<Scalar>::getAreas() const
195  {
196  return this->areas;
197  }
198 
199  template<typename Scalar>
200  void Form<Scalar>::setScalingFactor(double scalingFactor)
201  {
202  this->scaling_factor = scalingFactor;
203  }
204 
205  template<typename Scalar>
206  void Form<Scalar>::set_uExtOffset(int u_ext_offset)
207  {
208  this->u_ext_offset = u_ext_offset;
209  }
210 
211  template<typename Scalar>
212  void Form<Scalar>::set_ext(MeshFunction<Scalar>* ext)
213  {
214  this->ext.clear();
215  this->ext.push_back(ext);
216  }
217 
218  template<typename Scalar>
219  void Form<Scalar>::set_ext(Hermes::vector<MeshFunction<Scalar>*> ext)
220  {
221  this->ext = ext;
222  }
223 
224  template<typename Scalar>
225  Hermes::vector<MeshFunction<Scalar>*> Form<Scalar>::get_ext() const
226  {
227  return this->ext;
228  }
229 
230  template<typename Scalar>
231  MatrixForm<Scalar>::MatrixForm(unsigned int i, unsigned int j) :
232  Form<Scalar>(), sym(HERMES_NONSYM), i(i), j(j), previous_iteration_space_index(-1)
233  {
234  }
235 
236  template<typename Scalar>
237  MatrixForm<Scalar>::~MatrixForm()
238  {
239  }
240 
241  template<typename Scalar>
242  Scalar MatrixForm<Scalar>::value(int n, double *wt, Func<Scalar> *u_ext[], Func<double> *u, Func<double> *v,
243  Geom<double> *e, Func<Scalar> **ext) const
244  {
245  throw Hermes::Exceptions::MethodNotOverridenException("MatrixForm<Scalar>::value");
246  return 0.0;
247  }
248 
249  template<typename Scalar>
250  Hermes::Ord MatrixForm<Scalar>::ord(int n, double *wt, Func<Hermes::Ord> *u_ext[], Func<Hermes::Ord> *u, Func<Hermes::Ord> *v,
251  Geom<Hermes::Ord> *e, Func<Ord> **ext) const
252  {
253  throw Hermes::Exceptions::MethodNotOverridenException("MatrixForm<Scalar>::ord");
254  return Hermes::Ord();
255  }
256 
257  template<typename Scalar>
258  MatrixFormVol<Scalar>::MatrixFormVol(unsigned int i, unsigned int j) :
259  MatrixForm<Scalar>(i, j)
260  {
261  }
262 
263  template<typename Scalar>
264  MatrixFormVol<Scalar>::~MatrixFormVol()
265  {
266  }
267 
268  template<typename Scalar>
269  void MatrixFormVol<Scalar>::setSymFlag(SymFlag sym)
270  {
271  this->sym = sym;
272  }
273 
274  template<typename Scalar>
275  SymFlag MatrixFormVol<Scalar>::getSymFlag() const
276  {
277  return this->sym;
278  }
279 
280  template<typename Scalar>
281  MatrixFormVol<Scalar>* MatrixFormVol<Scalar>::clone() const
282  {
283  throw Hermes::Exceptions::MethodNotOverridenException("MatrixFormVol<Scalar>::clone()");
284  return NULL;
285  }
286 
287  template<typename Scalar>
288  MatrixFormSurf<Scalar>::MatrixFormSurf(unsigned int i, unsigned int j) :
289  MatrixForm<Scalar>(i, j)
290  {
291  }
292 
293  template<typename Scalar>
294  MatrixFormSurf<Scalar>* MatrixFormSurf<Scalar>::clone() const
295  {
296  throw Hermes::Exceptions::MethodNotOverridenException("MatrixFormSurf<Scalar>::clone()");
297  return NULL;
298  }
299 
300  template<typename Scalar>
301  MatrixFormSurf<Scalar>::~MatrixFormSurf()
302  {
303  }
304 
305  template<typename Scalar>
306  MatrixFormDG<Scalar>::MatrixFormDG(unsigned int i, unsigned int j) :
307  MatrixForm<Scalar>(i, j)
308  {
309  this->set_area(H2D_DG_INNER_EDGE);
310  }
311 
312  template<typename Scalar>
313  MatrixFormDG<Scalar>::~MatrixFormDG()
314  {
315  }
316 
317  template<typename Scalar>
318  MatrixFormDG<Scalar>* MatrixFormDG<Scalar>::clone() const
319  {
320  throw Hermes::Exceptions::MethodNotOverridenException("MatrixFormDG<Scalar>::clone()");
321  return NULL;
322  }
323 
324  template<typename Scalar>
325  VectorForm<Scalar>::VectorForm(unsigned int i) :
326  Form<Scalar>(), i(i)
327  {
328  }
329 
330  template<typename Scalar>
331  VectorForm<Scalar>::~VectorForm()
332  {
333  }
334 
335  template<typename Scalar>
336  VectorFormVol<Scalar>::VectorFormVol(unsigned int i) :
337  VectorForm<Scalar>(i)
338  {
339  }
340 
341  template<typename Scalar>
342  VectorFormVol<Scalar>::~VectorFormVol()
343  {
344  }
345 
346  template<typename Scalar>
347  Scalar VectorForm<Scalar>::value(int n, double *wt, Func<Scalar> *u_ext[], Func<double> *v,
348  Geom<double> *e, Func<Scalar> **ext) const
349  {
350  throw Hermes::Exceptions::MethodNotOverridenException("VectorForm<Scalar>::value");
351  return 0.0;
352  }
353 
354  template<typename Scalar>
355  Hermes::Ord VectorForm<Scalar>::ord(int n, double *wt, Func<Hermes::Ord> *u_ext[], Func<Hermes::Ord> *v,
356  Geom<Hermes::Ord> *e, Func<Ord> **ext) const
357  {
358  throw Hermes::Exceptions::MethodNotOverridenException("VectorForm<Scalar>::ord");
359  return Hermes::Ord();
360  }
361 
362  template<typename Scalar>
363  VectorFormVol<Scalar>* VectorFormVol<Scalar>::clone() const
364  {
365  throw Hermes::Exceptions::MethodNotOverridenException("VectorFormVol<Scalar>::clone()");
366  return NULL;
367  }
368 
369  template<typename Scalar>
370  VectorFormSurf<Scalar>::VectorFormSurf(unsigned int i) :
371  VectorForm<Scalar>(i)
372  {
373  }
374 
375  template<typename Scalar>
376  VectorFormSurf<Scalar>::~VectorFormSurf()
377  {
378  }
379 
380  template<typename Scalar>
381  VectorFormSurf<Scalar>* VectorFormSurf<Scalar>::clone() const
382  {
383  throw Hermes::Exceptions::MethodNotOverridenException("VectorFormSurf<Scalar>::clone()");
384  return NULL;
385  }
386 
387  template<typename Scalar>
388  VectorFormDG<Scalar>::VectorFormDG(unsigned int i) :
389  VectorForm<Scalar>(i)
390  {
391  this->set_area(H2D_DG_INNER_EDGE);
392  }
393 
394  template<typename Scalar>
395  VectorFormDG<Scalar>::~VectorFormDG()
396  {
397  }
398 
399  template<typename Scalar>
400  VectorFormDG<Scalar>* VectorFormDG<Scalar>::clone() const
401  {
402  throw Hermes::Exceptions::MethodNotOverridenException("VectorFormDG<Scalar>::clone()");
403  return NULL;
404  }
405 
406  template<typename Scalar>
407  void WeakForm<Scalar>::add_matrix_form(MatrixFormVol<Scalar>* form)
408  {
409  if(form->i >= neq || form->j >= neq)
410  throw Hermes::Exceptions::Exception("Invalid equation number.");
411  if(form->sym < -1 || form->sym > 1)
412  throw Hermes::Exceptions::Exception("\"sym\" must be -1, 0 or 1.");
413  if(form->sym < 0 && form->i == form->j)
414  throw Hermes::Exceptions::Exception("Only off-diagonal forms can be antisymmetric.");
415  if(mfvol.size() > 100)
416  {
417  this->warn("Large number of forms (> 100). Is this the intent?");
418  }
419 
420  form->set_weakform(this);
421  mfvol.push_back(form);
422  forms.push_back(form);
423  }
424 
425  template<typename Scalar>
426  void WeakForm<Scalar>::add_matrix_form_surf(MatrixFormSurf<Scalar>* form)
427  {
428  if(form->i >= neq || form->j >= neq)
429  throw Hermes::Exceptions::Exception("Invalid equation number.");
430 
431  form->set_weakform(this);
432  mfsurf.push_back(form);
433  forms.push_back(form);
434  }
435 
436  template<typename Scalar>
437  void WeakForm<Scalar>::add_matrix_form_DG(MatrixFormDG<Scalar>* form)
438  {
439  if(form->i >= neq || form->j >= neq)
440  throw Hermes::Exceptions::Exception("Invalid equation number.");
441 
442  form->set_weakform(this);
443  mfDG.push_back(form);
444  forms.push_back(form);
445  }
446 
447  template<typename Scalar>
448  void WeakForm<Scalar>::add_vector_form(VectorFormVol<Scalar>* form)
449  {
450  if(form->i >= neq)
451  throw Hermes::Exceptions::Exception("Invalid equation number.");
452  form->set_weakform(this);
453  vfvol.push_back(form);
454  forms.push_back(form);
455  }
456 
457  template<typename Scalar>
458  void WeakForm<Scalar>::add_vector_form_surf(VectorFormSurf<Scalar>* form)
459  {
460  if(form->i >= neq)
461  throw Hermes::Exceptions::Exception("Invalid equation number.");
462 
463  form->set_weakform(this);
464  vfsurf.push_back(form);
465  forms.push_back(form);
466  }
467 
468  template<typename Scalar>
469  void WeakForm<Scalar>::add_vector_form_DG(VectorFormDG<Scalar>* form)
470  {
471  if(form->i >= neq)
472  throw Hermes::Exceptions::Exception("Invalid equation number.");
473 
474  form->set_weakform(this);
475  vfDG.push_back(form);
476  forms.push_back(form);
477  }
478 
479  template<typename Scalar>
480  Hermes::vector<Form<Scalar> *> WeakForm<Scalar>::get_forms() const
481  {
482  return forms;
483  }
484 
485  template<typename Scalar>
486  Hermes::vector<MatrixFormVol<Scalar> *> WeakForm<Scalar>::get_mfvol() const
487  {
488  return mfvol;
489  }
490  template<typename Scalar>
491  Hermes::vector<MatrixFormSurf<Scalar> *> WeakForm<Scalar>::get_mfsurf() const
492  {
493  return mfsurf;
494  }
495  template<typename Scalar>
496  Hermes::vector<MatrixFormDG<Scalar> *> WeakForm<Scalar>::get_mfDG() const
497  {
498  return mfDG;
499  }
500  template<typename Scalar>
501  Hermes::vector<VectorFormVol<Scalar> *> WeakForm<Scalar>::get_vfvol() const
502  {
503  return vfvol;
504  }
505  template<typename Scalar>
506  Hermes::vector<VectorFormSurf<Scalar> *> WeakForm<Scalar>::get_vfsurf() const
507  {
508  return vfsurf;
509  }
510  template<typename Scalar>
511  Hermes::vector<VectorFormDG<Scalar> *> WeakForm<Scalar>::get_vfDG() const
512  {
513  return vfDG;
514  }
515 
516  template<typename Scalar>
517  bool** WeakForm<Scalar>::get_blocks(bool force_diagonal_blocks) const
518  {
519  bool** blocks = new_matrix<bool>(neq, neq);
520  for (unsigned int i = 0; i < neq; i++)
521  {
522  for (unsigned int j = 0; j < neq; j++)
523  blocks[i][j] = false;
524  if(force_diagonal_blocks)
525  blocks[i][i] = true;
526  }
527  for (unsigned i = 0; i < mfvol.size(); i++)
528  {
529  if(fabs(mfvol[i]->scaling_factor) > 1e-12)
530  blocks[mfvol[i]->i][mfvol[i]->j] = true;
531  if(mfvol[i]->sym)
532  if(fabs(mfvol[i]->scaling_factor) > 1e-12)
533  blocks[mfvol[i]->j][mfvol[i]->i] = true;
534  }
535  for (unsigned i = 0; i < mfsurf.size(); i++)
536  {
537  if(fabs(mfsurf[i]->scaling_factor) > 1e-12)
538  blocks[mfsurf[i]->i][mfsurf[i]->j] = true;
539  }
540 
541  return blocks;
542  }
543 
544  template<typename Scalar>
545  void WeakForm<Scalar>::set_current_time(double time)
546  {
547  current_time = time;
548  }
549 
550  template<typename Scalar>
551  double WeakForm<Scalar>::get_current_time() const
552  {
553  return current_time;
554  }
555 
556  template<typename Scalar>
557  void WeakForm<Scalar>::set_current_time_step(double time_step)
558  {
559  current_time_step = time_step;
560  }
561 
562  template<typename Scalar>
563  double WeakForm<Scalar>::get_current_time_step() const
564  {
565  return current_time_step;
566  }
567 
568  template class HERMES_API WeakForm<double>;
569  template class HERMES_API WeakForm<std::complex<double> >;
570  template class HERMES_API Form<double>;
571  template class HERMES_API Form<std::complex<double> >;
572  template class HERMES_API MatrixForm<double>;
573  template class HERMES_API MatrixForm<std::complex<double> >;
574  template class HERMES_API MatrixFormVol<double>;
575  template class HERMES_API MatrixFormVol<std::complex<double> >;
576  template class HERMES_API MatrixFormSurf<double>;
577  template class HERMES_API MatrixFormSurf<std::complex<double> >;
578  template class HERMES_API MatrixFormDG<double>;
579  template class HERMES_API MatrixFormDG<std::complex<double> >;
580  template class HERMES_API VectorForm<double>;
581  template class HERMES_API VectorForm<std::complex<double> >;
582  template class HERMES_API VectorFormVol<double>;
583  template class HERMES_API VectorFormVol<std::complex<double> >;
584  template class HERMES_API VectorFormSurf<double>;
585  template class HERMES_API VectorFormSurf<std::complex<double> >;
586  template class HERMES_API VectorFormDG<double>;
587  template class HERMES_API VectorFormDG<std::complex<double> >;
588  }
589 }