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"
33 template<
typename Scalar>
38 template<
typename Scalar>
39 WeakFormSharedPtr<Scalar>::WeakFormSharedPtr(
const WeakFormSharedPtr& other) :
std::tr1::shared_ptr<
Hermes::Hermes2D::WeakForm<Scalar> >(other)
43 template<
typename Scalar>
44 void WeakFormSharedPtr<Scalar>::operator=(
const WeakFormSharedPtr& other)
46 std::tr1::shared_ptr<Hermes::Hermes2D::WeakForm<Scalar> >::operator=(other);
49 static bool warned_nonOverride =
false;
54 static const std::string H2D_DG_INNER_EDGE =
"-1234567";
56 template<
typename Scalar>
61 this->is_matfree = mat_free;
64 template<
typename Scalar>
70 template<
typename Scalar>
73 for (
unsigned int i = 0; i < this->forms.size(); i++)
74 delete get_forms()[i];
78 template<
typename Scalar>
81 if (!warned_nonOverride)
82 #pragma omp critical (warning_weakform_nonOverride)
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;
91 template<
typename Scalar>
97 template<
typename Scalar>
103 template<
typename Scalar>
109 template<
typename Scalar>
113 this->mfsurf.clear();
116 this->vfsurf.clear();
121 for (
unsigned int i = 0; i < other_wf->forms.size(); i++)
126 form->copy_base(other_wf->forms[i]);
127 this->forms.push_back(form);
131 if (
dynamic_cast<VectorFormVol<Scalar>*
>(other_wf->forms[i]) !=
nullptr)
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()));
139 if (
dynamic_cast<MatrixFormSurf<Scalar>*
>(other_wf->forms[i]) !=
nullptr)
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()));
147 if (
dynamic_cast<VectorFormSurf<Scalar>*
>(other_wf->forms[i]) !=
nullptr)
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()));
155 if (
dynamic_cast<MatrixFormDG<Scalar>*
>(other_wf->forms[i]) !=
nullptr)
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()));
163 if (
dynamic_cast<VectorFormDG<Scalar>*
>(other_wf->forms[i]) !=
nullptr)
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()));
173 for (
unsigned int i = 0; i < other_wf->forms.size(); i++)
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;
179 this->cloneMemberExtFunctions(other_wf->ext, this->ext);
180 this->u_ext_fn = other_wf->u_ext_fn;
183 template<
typename Scalar>
184 void WeakForm<Scalar>::cloneMemberExtFunctions(std::vector<MeshFunctionSharedPtr<Scalar> > source_ext, std::vector<MeshFunctionSharedPtr<Scalar> >& cloned_ext)
187 for (
unsigned int i = 0; i < source_ext.size(); i++)
189 Solution<Scalar>* originalSln =
dynamic_cast<Solution<Scalar>*
>(source_ext[i].get());
192 Solution<Scalar>* newSln =
nullptr;
193 if (originalSln->get_type() == HERMES_SLN)
195 newSln =
new Solution < Scalar > ;
196 newSln->copy(source_ext[i].
get());
199 newSln =
static_cast<Solution<Scalar>*
>(originalSln->clone());
201 newSln->set_type(originalSln->get_type());
202 cloned_ext.push_back(newSln);
205 cloned_ext.push_back(source_ext[i]->clone());
209 template<
typename Scalar>
221 template<
typename Scalar>
225 this->ext.push_back(ext);
228 template<
typename Scalar>
234 template<
typename Scalar>
237 this->u_ext_fn.clear();
238 this->u_ext_fn.push_back(ext);
241 template<
typename Scalar>
244 this->u_ext_fn = ext;
247 template<
typename Scalar>
248 template<
typename FormType>
251 for (
unsigned short form_i = 0; form_i < forms_to_process.size(); form_i++)
255 for (
unsigned short marker_i = 0; marker_i < form->
areas.size(); marker_i++)
257 if (form->
areas[marker_i] == HERMES_ANY)
264 Mesh::MarkersConversion::IntValid marker;
266 marker = spaces[form->i]->get_mesh()->get_boundary_markers_conversion().get_internal_marker(form->
areas[marker_i]);
268 marker = spaces[form->i]->get_mesh()->get_element_markers_conversion().get_internal_marker(form->
areas[marker_i]);
273 throw Exceptions::Exception(
"Marker not valid in assembling: %s.", form->
areas[marker_i].c_str());
278 template<
typename Scalar>
279 void WeakForm<Scalar>::processFormMarkers(
const std::vector<SpaceSharedPtr<Scalar> > spaces)
281 processFormMarkers(spaces,
false, this->mfvol);
282 processFormMarkers(spaces,
false, this->vfvol);
283 processFormMarkers(spaces,
true, this->mfsurf);
284 processFormMarkers(spaces,
true, this->vfsurf);
287 template<
typename Scalar>
288 bool WeakForm<Scalar>::is_DG()
const
290 if (this->mfDG.empty() && this->vfDG.empty())
295 template<
typename Scalar>
301 template<
typename Scalar>
304 areas.push_back(HERMES_ANY);
308 template<
typename Scalar>
313 template<
typename Scalar>
314 void Form<Scalar>::set_current_stage_time(
double time)
319 template<
typename Scalar>
320 double Form<Scalar>::get_current_stage_time()
const
325 template<
typename Scalar>
329 areas.push_back(area);
331 template<
typename Scalar>
337 template<
typename Scalar>
338 std::vector<std::string> Form<Scalar>::getAreas()
const
343 template<
typename Scalar>
346 this->scaling_factor = scalingFactor;
349 template<
typename Scalar>
353 this->ext.push_back(ext);
356 template<
typename Scalar>
362 template<
typename Scalar>
365 this->u_ext_fn.clear();
366 this->u_ext_fn.push_back(ext);
369 template<
typename Scalar>
370 void Form<Scalar>::set_u_ext_fn(std::vector<UExtFunctionSharedPtr<Scalar> > ext)
372 this->u_ext_fn = ext;
375 template<
typename Scalar>
376 std::vector<MeshFunctionSharedPtr<Scalar> > Form<Scalar>::get_ext()
const
381 template<
typename Scalar>
382 void Form<Scalar>::copy_base(Form<Scalar>* other_form)
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;
390 template<
typename Scalar>
392 Form<Scalar>(i), sym(HERMES_NONSYM), j(j)
397 template<
typename Scalar>
402 template<
typename Scalar>
408 template<
typename Scalar>
413 template<
typename Scalar>
414 void MatrixFormVol<Scalar>::setSymFlag(
SymFlag sym)
419 template<
typename Scalar>
420 SymFlag MatrixFormVol<Scalar>::getSymFlag()
const
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
429 throw Hermes::Exceptions::MethodNotOverridenException(
"MatrixFormVol<Scalar>::value");
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
437 throw Hermes::Exceptions::MethodNotOverridenException(
"MatrixFormVol<Scalar>::ord");
438 return Hermes::Ord();
441 template<
typename Scalar>
442 MatrixFormVol<Scalar>* MatrixFormVol<Scalar>::clone()
const
444 throw Hermes::Exceptions::MethodNotOverridenException(
"MatrixFormVol<Scalar>::clone()");
448 template<
typename Scalar>
454 template<
typename Scalar>
458 throw Hermes::Exceptions::MethodNotOverridenException(
"MatrixFormSurf<Scalar>::value");
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
466 throw Hermes::Exceptions::MethodNotOverridenException(
"MatrixFormSurf<Scalar>::ord");
467 return Hermes::Ord();
470 template<
typename Scalar>
471 MatrixFormSurf<Scalar>* MatrixFormSurf<Scalar>::clone()
const
473 throw Hermes::Exceptions::MethodNotOverridenException(
"MatrixFormSurf<Scalar>::clone()");
477 template<
typename Scalar>
478 MatrixFormSurf<Scalar>::~MatrixFormSurf()
482 template<
typename Scalar>
484 Form<Scalar>(i), j(j)
490 template<
typename Scalar>
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
499 throw Hermes::Exceptions::MethodNotOverridenException(
"MatrixFormDG<Scalar>::value");
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
507 throw Hermes::Exceptions::MethodNotOverridenException(
"MatrixFormDG<Scalar>::ord");
508 return Hermes::Ord();
511 template<
typename Scalar>
512 MatrixFormDG<Scalar>* MatrixFormDG<Scalar>::clone()
const
514 throw Hermes::Exceptions::MethodNotOverridenException(
"MatrixFormDG<Scalar>::clone()");
518 template<
typename Scalar>
525 template<
typename Scalar>
530 template<
typename Scalar>
536 template<
typename Scalar>
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
545 throw Hermes::Exceptions::MethodNotOverridenException(
"VectorFormVol<Scalar>::value");
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
553 throw Hermes::Exceptions::MethodNotOverridenException(
"VectorFormVol<Scalar>::ord");
554 return Hermes::Ord();
557 template<
typename Scalar>
558 VectorFormVol<Scalar>* VectorFormVol<Scalar>::clone()
const
560 throw Hermes::Exceptions::MethodNotOverridenException(
"VectorFormVol<Scalar>::clone()");
564 template<
typename Scalar>
570 template<
typename Scalar>
574 throw Hermes::Exceptions::MethodNotOverridenException(
"VectorFormSurf<Scalar>::value");
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
582 throw Hermes::Exceptions::MethodNotOverridenException(
"VectorFormSurf<Scalar>::ord");
583 return Hermes::Ord();
586 template<
typename Scalar>
587 VectorFormSurf<Scalar>::~VectorFormSurf()
591 template<
typename Scalar>
592 VectorFormSurf<Scalar>* VectorFormSurf<Scalar>::clone()
const
594 throw Hermes::Exceptions::MethodNotOverridenException(
"VectorFormSurf<Scalar>::clone()");
598 template<
typename Scalar>
605 template<
typename Scalar>
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
614 throw Hermes::Exceptions::MethodNotOverridenException(
"VectorFormDG<Scalar>::value");
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
622 throw Hermes::Exceptions::MethodNotOverridenException(
"VectorFormDG<Scalar>::ord");
623 return Hermes::Ord();
626 template<
typename Scalar>
627 VectorFormDG<Scalar>* VectorFormDG<Scalar>::clone()
const
629 throw Hermes::Exceptions::MethodNotOverridenException(
"VectorFormDG<Scalar>::clone()");
633 template<
typename Scalar>
638 this->previous_iteration_space_index %= wf->
original_neq;
641 template<
typename Scalar>
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)
652 this->warn(
"Large number of forms (> 100). Is this the intent?");
656 mfvol.push_back(form);
657 forms.push_back(form);
660 template<
typename Scalar>
663 if (form->i >= neq || form->j >= neq)
664 throw Hermes::Exceptions::Exception(
"Invalid equation number.");
667 mfsurf.push_back(form);
668 forms.push_back(form);
671 template<
typename Scalar>
674 if (form->i >= neq || form->j >= neq)
675 throw Hermes::Exceptions::Exception(
"Invalid equation number.");
678 mfDG.push_back(form);
679 forms.push_back(form);
682 template<
typename Scalar>
686 throw Hermes::Exceptions::Exception(
"Invalid equation number.");
688 vfvol.push_back(form);
689 forms.push_back(form);
692 template<
typename Scalar>
696 throw Hermes::Exceptions::Exception(
"Invalid equation number.");
699 vfsurf.push_back(form);
700 forms.push_back(form);
703 template<
typename Scalar>
707 throw Hermes::Exceptions::Exception(
"Invalid equation number.");
710 vfDG.push_back(form);
711 forms.push_back(form);
714 template<
typename Scalar>
720 template<
typename Scalar>
725 template<
typename Scalar>
726 std::vector<MatrixFormSurf<Scalar> *> WeakForm<Scalar>::get_mfsurf()
const
730 template<
typename Scalar>
731 std::vector<MatrixFormDG<Scalar> *> WeakForm<Scalar>::get_mfDG()
const
735 template<
typename Scalar>
736 std::vector<VectorFormVol<Scalar> *> WeakForm<Scalar>::get_vfvol()
const
740 template<
typename Scalar>
741 std::vector<VectorFormSurf<Scalar> *> WeakForm<Scalar>::get_vfsurf()
const
745 template<
typename Scalar>
746 std::vector<VectorFormDG<Scalar> *> WeakForm<Scalar>::get_vfDG()
const
751 template<
typename Scalar>
752 bool** WeakForm<Scalar>::get_blocks(
bool force_diagonal_blocks)
const
754 bool** blocks = new_matrix<bool>(neq, neq);
755 for (
unsigned int i = 0; i < neq; i++)
757 for (
unsigned int j = 0; j < neq; j++)
758 blocks[i][j] =
false;
759 if (force_diagonal_blocks)
762 for (
unsigned i = 0; i < mfvol.size(); i++)
764 if (fabs(mfvol[i]->scaling_factor) > Hermes::HermesSqrtEpsilon)
765 blocks[mfvol[i]->i][mfvol[i]->j] =
true;
767 if (fabs(mfvol[i]->scaling_factor) > Hermes::HermesSqrtEpsilon)
768 blocks[mfvol[i]->j][mfvol[i]->i] =
true;
770 for (
unsigned i = 0; i < mfsurf.size(); i++)
772 if (fabs(mfsurf[i]->scaling_factor) > Hermes::HermesSqrtEpsilon)
773 blocks[mfsurf[i]->i][mfsurf[i]->j] =
true;
775 for (
unsigned i = 0; i < mfDG.size(); i++)
777 if (fabs(mfDG[i]->scaling_factor) > Hermes::HermesSqrtEpsilon)
778 blocks[mfDG[i]->i][mfDG[i]->j] =
true;
784 template<
typename Scalar>
790 template<
typename Scalar>
796 template<
typename Scalar>
799 current_time_step = time_step;
802 template<
typename Scalar>
805 return current_time_step;
Stores one element of a mesh.
Used to pass the instances of Space around.
Common definitions for Hermes2D.
::xsd::cxx::tree::time< char, simple_type > time
C++ type corresponding to the time XML Schema built-in type.
SymFlag
Bilinear form symmetry flag, see WeakForm::add_matrix_form.
Calculated function values (from the class Function) on an element for assembling.
::xsd::cxx::tree::string< char, simple_type > string
C++ type corresponding to the string XML Schema built-in type.