16 #include "weakform_library/weakforms_neutronics.h"
17 #include "weakform_library/integrals_h1.h"
25 namespace WeakFormsNeutronics
27 namespace Monoenergetic
31 template<
typename Scalar>
32 DefaultWeakFormFixedSource<Scalar>::DefaultWeakFormFixedSource(std::vector<std::string> regions,
33 std::vector<double> D_map,
34 std::vector<double> Sigma_a_map,
35 std::vector<double> Q_map) : WeakForm<Scalar>(1)
37 using namespace WeakFormsH1;
39 for (
unsigned int i = 0; i < regions.size(); i++)
43 this->add_matrix_form(
new DefaultJacobianDiffusion<Scalar>(0, 0, regions[i],
new Hermes1DFunction<Scalar>(D_map[i]),
46 this->add_matrix_form(
new DefaultMatrixFormVol<Scalar>(0, 0, regions[i],
new Hermes2DFunction<Scalar>(Sigma_a_map[i]),
51 this->add_vector_form(
new DefaultResidualDiffusion<Scalar>(0, regions[i],
new Hermes1DFunction<Scalar>(D_map[i])));
53 this->add_vector_form(
new DefaultResidualVol<Scalar>(0, regions[i],
new Hermes2DFunction<Scalar>(Sigma_a_map[i])));
55 this->add_vector_form(
new DefaultVectorFormVol<Scalar>(0, regions[i],
new Hermes2DFunction<Scalar>(-Q_map[i])));
63 namespace MaterialProperties
67 void MaterialPropertyMaps::extend_to_multigroup(
const MaterialPropertyMap0& mrsg_map,
68 MaterialPropertyMap1 *mrmg_map)
71 this->warn(W_MG_EXTENSION);
73 MaterialPropertyMap0::const_iterator it;
74 for (it = mrsg_map.begin(); it != mrsg_map.end(); ++it)
75 (*mrmg_map)[it->first].assign(G, it->second);
78 void MaterialPropertyMaps::extend_to_multiregion(
const rank1& srmg_array,
79 MaterialPropertyMap1 *mrmg_map)
81 if (materials_list.empty())
82 throw Hermes::Exceptions::Exception(E_MR_EXTENSION);
84 std::set<std::string>::const_iterator it;
85 for (it = materials_list.begin(); it != materials_list.end(); ++it)
86 (*mrmg_map)[*it] = srmg_array;
89 void MaterialPropertyMaps::extend_to_multiregion_multigroup(
const rank0& srsg_value,
90 MaterialPropertyMap1 *mrmg_map)
92 if (materials_list.empty())
93 throw Hermes::Exceptions::Exception(E_MR_EXTENSION);
95 this->warn(W_MG_EXTENSION);
97 std::set<std::string>::const_iterator it;
98 for (it = materials_list.begin(); it != materials_list.end(); ++it)
99 (*mrmg_map)[*it].assign(G, srsg_value);
102 void MaterialPropertyMaps::fill_with(
double c, MaterialPropertyMap1 *mrmg_map)
104 if (materials_list.empty())
105 throw Hermes::Exceptions::Exception(E_MR_EXTENSION);
107 std::set<std::string>::const_iterator it;
108 for (it = materials_list.begin(); it != materials_list.end(); ++it)
109 (*mrmg_map)[*it].assign(G, c);
112 void MaterialPropertyMaps::validate()
114 using namespace ValidationFunctors;
116 if (fission_multigroup_structure.empty())
117 fission_multigroup_structure = bool1(G,
true);
121 fill_with(0.0, &chi);
122 MaterialPropertyMap1::iterator it = chi.begin();
123 for (; it != chi.end(); ++it)
125 fission_multigroup_structure = bool1(G,
false);
126 fission_multigroup_structure[0] =
true;
129 if (nu.empty() && !nuSigma_f.empty() && !Sigma_f.empty())
130 nu = NDArrayMapOp::divide<rank1>(nuSigma_f, Sigma_f);
131 else if (nuSigma_f.empty() && !nu.empty() && !Sigma_f.empty())
132 nuSigma_f = NDArrayMapOp::multiply<rank1>(nu, Sigma_f);
133 else if (Sigma_f.empty() && !nuSigma_f.empty() && !nu.empty())
134 Sigma_f = NDArrayMapOp::divide<rank1>(nuSigma_f, nu);
135 else if (!Sigma_f.empty() && !nuSigma_f.empty() && !nu.empty())
137 MaterialPropertyMap1 diff = NDArrayMapOp::subtract<rank1>(nuSigma_f,
138 NDArrayMapOp::multiply<rank1>(nu, Sigma_f));
139 std::for_each(diff.begin(), diff.end(), ensure_trivial());
143 this->warn(W_NO_FISSION);
145 fill_with(0.0, &chi);
146 fill_with(0.0, &Sigma_f);
149 if ((nu.size() != Sigma_f.size()) || (nu.size() != chi.size()))
150 throw Hermes::Exceptions::Exception(E_NONMATCHING_PROPERTIES);
152 if (Sigma_f.size() > 0)
154 std::for_each(nu.begin(), nu.end(), ensure_size(G));
155 std::for_each(Sigma_f.begin(), Sigma_f.end(), ensure_size(G));
156 std::for_each(chi.begin(), chi.end(), ensure_size(G));
159 if (Sigma_a.size() > 0)
164 MaterialPropertyMap1::const_iterator ita = Sigma_a.begin();
165 MaterialPropertyMap1::const_iterator itf = Sigma_f.begin();
166 for (; ita != Sigma_a.end(); ++ita, ++itf)
168 rank1::const_iterator a = ita->second.begin();
169 rank1::const_iterator f = itf->second.begin();
171 for (; a != ita->second.end(); ++a, ++f)
173 this->warn(W_SA_LT_SF);
178 const rank1& MaterialPropertyMaps::get_Sigma_f(
std::string material)
const
180 if (material ==
"-999")
return this->Sigma_f.begin()->second;
184 MaterialPropertyMap1::const_iterator data = this->Sigma_f.find(material);
185 if (data != this->Sigma_f.end())
189 throw Hermes::Exceptions::Exception(E_INVALID_MARKER);
191 return *(
new rank1());
194 const rank1& MaterialPropertyMaps::get_nu(
std::string material)
const
196 if (material ==
"-999")
return this->nu.begin()->second;
198 MaterialPropertyMap1::const_iterator data = this->nu.find(material);
199 if (data != this->nu.end())
203 throw Hermes::Exceptions::Exception(E_INVALID_MARKER);
205 return *(
new rank1());
208 const rank1& MaterialPropertyMaps::get_chi(
std::string material)
const
210 if (material ==
"-999")
return this->chi.begin()->second;
212 MaterialPropertyMap1::const_iterator data = this->chi.find(material);
213 if (data != this->chi.end())
217 throw Hermes::Exceptions::Exception(E_INVALID_MARKER);
219 return *(
new rank1());
223 std::ostream & operator<< (std::ostream& os,
const MaterialPropertyMaps& matprop)
228 os << setw(12) <<
"target group" << setw(10) <<
"chi" << setw(10) <<
"nu";
229 os << setw(10) <<
"Sigma_f" << std::endl;
231 MaterialPropertyMap1::const_iterator data_elem = matprop.chi.begin();
232 for (; data_elem != matprop.chi.end(); ++data_elem)
234 string mat = data_elem->first;
236 os << setw(80) << setfill(
'-') <<
' ' << std::endl << setfill(
' ');
237 os << setw(40) << mat << std::endl;
238 os << setw(80) << setfill(
'-') <<
' ' << std::endl << setfill(
' ');
239 for (
unsigned int gto = 0; gto < matprop.G; gto++)
241 os << setw(6) << gto << setw(6) <<
' ';
242 os << setw(10) << matprop.get_chi(mat)[gto];
243 os << setw(10) << matprop.get_nu(mat)[gto];
244 os << setw(10) << matprop.get_Sigma_f(mat)[gto];
257 MaterialPropertyMap1 MaterialPropertyMaps::extract_map2_diagonals(
const MaterialPropertyMap2& map2)
259 MaterialPropertyMap1 diags;
261 MaterialPropertyMap2::const_iterator map2_it = map2.begin();
262 for (; map2_it != map2.end(); ++map2_it)
264 diags[map2_it->first].reserve(G);
265 for (
unsigned int g = 0; g < G; g++)
266 diags[map2_it->first].push_back(map2_it->second[g][g]);
272 MaterialPropertyMap1 MaterialPropertyMaps::sum_map2_columns(
const MaterialPropertyMap2& map2)
274 MaterialPropertyMap1 summed;
276 MaterialPropertyMap2::const_iterator map2_it = map2.begin();
277 for (; map2_it != map2.end(); ++map2_it)
279 summed[map2_it->first].reserve(G);
280 for (
unsigned int gfrom = 0; gfrom < G; gfrom++)
284 for (
unsigned int gto = 0; gto < G; gto++)
285 sum += map2_it->second[gto][gfrom];
287 summed[map2_it->first].push_back(sum);
294 MaterialPropertyMap2 MaterialPropertyMaps::create_map2_by_diagonals(
const MaterialPropertyMap1& diags)
296 MaterialPropertyMap2 map2;
298 MaterialPropertyMap1::const_iterator diags_it = diags.begin();
299 for (; diags_it != diags.end(); ++diags_it)
301 map2[diags_it->first].resize(G, rank1(G, 0.0));
303 for (
unsigned int g = 0; g < G; g++)
304 map2[diags_it->first][g][g] = diags_it->second[g];
310 void MaterialPropertyMaps::fill_with(
double c, MaterialPropertyMap2 *mrmg_map)
312 std::set<std::string>::const_iterator it;
313 for (it = materials_list.begin(); it != materials_list.end(); ++it)
314 (*mrmg_map)[*it].assign(G, rank1(G, c));
317 void MaterialPropertyMaps::validate()
319 Common::MaterialPropertyMaps::validate();
321 bool D_given = !D.empty();
322 bool Sigma_r_given = !Sigma_r.empty();
323 bool Sigma_s_given = !Sigma_s.empty();
324 bool Sigma_t_given = !Sigma_t.empty();
325 bool Sigma_a_given = !Sigma_a.empty();
326 bool Sigma_f_given = !Sigma_f.empty();
327 bool src_given = !src.empty();
340 Sigma_s = create_map2_by_diagonals(Common::NDArrayMapOp::subtract<rank1>(Sigma_t, Sigma_a));
346 this->warn(W_NO_SCATTERING);
347 fill_with(0.0, &Sigma_s);
350 Sigma_s_given =
true;
359 this->warn(W_NO_SCATTERING);
360 fill_with(0.0, &Sigma_s);
361 Sigma_s_given =
true;
365 Sigma_t = Common::NDArrayMapOp::add<rank1>(Sigma_a, sum_map2_columns(Sigma_s));
372 throw Hermes::Exceptions::Exception(E_INSUFFICIENT_DATA);
375 Sigma_t_given =
true;
378 Sigma_r = Common::NDArrayMapOp::subtract<rank1>(Sigma_t, extract_map2_diagonals(Sigma_s));
379 Sigma_r_given =
true;
384 if (scattering_multigroup_structure.empty())
385 scattering_multigroup_structure = bool2(G, std::vector<bool>(G,
true));
395 Sigma_s = create_map2_by_diagonals(Common::NDArrayMapOp::subtract<rank1>(Sigma_t, Sigma_r));
397 scattering_multigroup_structure = bool2(G, std::vector<bool>(G,
false));
398 for (
unsigned int gto = 0; gto < G; gto++)
399 for (
unsigned int gfrom = 0; gfrom < G; gfrom++)
401 scattering_multigroup_structure[gto][gfrom] =
true;
405 this->warn(W_NO_SCATTERING);
406 fill_with(0.0, &Sigma_s);
407 scattering_multigroup_structure = bool2(G, std::vector<bool>(G,
false));
410 Sigma_s_given =
true;
417 MaterialPropertyMap1::const_iterator Sr_elem = Sigma_r.begin();
418 for (; Sr_elem != Sigma_r.end(); ++Sr_elem)
419 for (
unsigned int g = 0; g < G; g++)
420 D[Sr_elem->first][g] = 1. / (3.*Sr_elem->second[g]);
425 if ((D.size() != Sigma_r.size()) || (D.size() != Sigma_s.size()) || (src_given && D.size() != src.size()))
426 throw Hermes::Exceptions::Exception(E_NONMATCHING_PROPERTIES);
428 using ValidationFunctors::ensure_size;
429 std::for_each(Sigma_s.begin(), Sigma_s.end(), ensure_size(G, G));
430 std::for_each(Sigma_r.begin(), Sigma_r.end(), ensure_size(G));
431 std::for_each(src.begin(), src.end(), ensure_size(G));
432 std::for_each(D.begin(), D.end(), ensure_size(G));
435 const rank2& MaterialPropertyMaps::get_Sigma_s(
std::string material)
const
437 if (material ==
"-999")
return this->Sigma_s.begin()->second;
441 MaterialPropertyMap2::const_iterator data = this->Sigma_s.find(material);
442 if (data != this->Sigma_s.end())
446 throw Hermes::Exceptions::Exception(E_INVALID_MARKER);
448 return *(
new rank2());
451 const rank1& MaterialPropertyMaps::get_Sigma_r(
std::string material)
const
453 if (material ==
"-999")
return this->Sigma_r.begin()->second;
455 MaterialPropertyMap1::const_iterator data = this->Sigma_r.find(material);
456 if (data != this->Sigma_r.end())
460 throw Hermes::Exceptions::Exception(E_INVALID_MARKER);
462 return *(
new rank1());
465 const rank1& MaterialPropertyMaps::get_D(
std::string material)
const
467 if (material ==
"-999")
return this->D.begin()->second;
469 MaterialPropertyMap1::const_iterator data = this->D.find(material);
470 if (data != this->D.end())
474 throw Hermes::Exceptions::Exception(E_INVALID_MARKER);
476 return *(
new rank1());
479 const rank1& MaterialPropertyMaps::get_src(
std::string material)
const
481 if (material ==
"-999")
return this->src.begin()->second;
483 MaterialPropertyMap1::const_iterator data = this->src.find(material);
484 if (data != this->src.end())
488 throw Hermes::Exceptions::Exception(E_INVALID_MARKER);
490 return *(
new rank1());
494 HERMES_API std::ostream & operator<< (std::ostream& os,
const MaterialPropertyMaps& matprop)
498 os << static_cast<const Common::MaterialPropertyMaps&>(matprop) << std::endl;
500 os << setw(12) <<
"target group" << setw(10) <<
"D" << setw(10) <<
"Sigma_r";
501 os << setw(10) <<
"ext. src" << setw(22) <<
"Sigma_s" << std::endl;
503 MaterialPropertyMap1::const_iterator data_elem = matprop.Sigma_r.begin();
504 for (; data_elem != matprop.Sigma_r.end(); ++data_elem)
506 string mat = data_elem->first;
508 os << setw(80) << setfill(
'-') <<
' ' << std::endl << setfill(
' ');
509 os << setw(40) << mat << std::endl;
510 os << setw(80) << setfill(
'-') <<
' ' << std::endl << setfill(
' ');
511 for (
unsigned int gto = 0; gto < matprop.G; gto++)
513 os << setw(6) << gto << setw(6) <<
' ';
514 os << setw(10) << matprop.get_D(mat)[gto];
515 os << setw(10) << matprop.get_Sigma_r(mat)[gto];
517 if (matprop.src.empty())
520 os << matprop.get_src(mat)[gto];
522 for (
unsigned int gfrom = 0; gfrom < matprop.G; gfrom++)
523 os << setw(8) << matprop.get_Sigma_s(mat)[gto][gfrom];
529 return os << std::endl;
534 namespace ElementaryForms
538 template<
typename ScalarClass>
539 template<
typename Real,
typename Scalar>
540 Scalar VacuumBoundaryCondition::Jacobian<ScalarClass>::matrix_form(
int n,
double *wt, Func<Scalar> *u_ext[], Func<Real> *u,
541 Func<Real> *v, GeomSurf<Real> *e, Func<Scalar> **ext)
const
545 if (geom_type == HERMES_PLANAR)
546 result = 0.5 * int_u_v<Real, Scalar>(n, wt, u, v);
547 else if (geom_type == HERMES_AXISYM_X)
548 result = 0.5 * int_y_u_v<Real, Scalar>(n, wt, u, v, e);
550 result = 0.5 * int_x_u_v<Real, Scalar>(n, wt, u, v, e);
555 template<
typename ScalarClass>
556 template<
typename Real,
typename Scalar>
557 Scalar VacuumBoundaryCondition::Residual<ScalarClass>::vector_form(
int n,
double *wt, Func<Scalar> *u_ext[],
558 Func<Real> *v, GeomSurf<Real> *e, Func<Scalar> **ext)
const
562 if (geom_type == HERMES_PLANAR)
563 result = 0.5 * int_u_ext_v<Real, Scalar>(n, wt, u_ext[g], v);
564 else if (geom_type == HERMES_AXISYM_X)
565 result = 0.5 * int_y_u_ext_v<Real, Scalar>(n, wt, u_ext[g], v, e);
567 result = 0.5 * int_x_u_ext_v<Real, Scalar>(n, wt, u_ext[g], v, e);
572 template<
typename ScalarClass>
573 template<
typename Real,
typename Scalar>
574 Scalar DiffusionReaction::Jacobian<ScalarClass>::matrix_form(
int n,
double *wt, Func<Scalar> *u_ext[], Func<Real> *u,
575 Func<Real> *v, GeomVol<Real> *e, Func<Scalar> **ext)
const
579 std::string mat = this->mesh->get_element_markers_conversion().get_user_marker(e->elem_marker).marker;
580 rank1 D_elem = matprop.get_D(mat);
581 rank1 Sigma_r_elem = matprop.get_Sigma_r(mat);
583 if (geom_type == HERMES_PLANAR)
585 result = D_elem[g] * int_grad_u_grad_v<Real, Scalar>(n, wt, u, v) +
586 Sigma_r_elem[g] * int_u_v<Real, Scalar>(n, wt, u, v);
590 if (geom_type == HERMES_AXISYM_X)
592 result = D_elem[g] * int_y_grad_u_grad_v<Real, Scalar>(n, wt, u, v, e) +
593 Sigma_r_elem[g] * int_y_u_v<Real, Scalar>(n, wt, u, v, e);
597 result = D_elem[g] * int_x_grad_u_grad_v<Real, Scalar>(n, wt, u, v, e) +
598 Sigma_r_elem[g] * int_x_u_v<Real, Scalar>(n, wt, u, v, e);
604 template<
typename ScalarClass>
605 template<
typename Real,
typename Scalar>
606 Scalar DiffusionReaction::Residual<ScalarClass>::vector_form(
int n,
double *wt, Func<Scalar> *u_ext[],
607 Func<Real> *v, GeomVol<Real> *e, Func<Scalar> **ext)
const
611 std::string mat = this->mesh->get_element_markers_conversion().get_user_marker(e->elem_marker).marker;
612 rank1 D_elem = matprop.get_D(mat);
613 rank1 Sigma_r_elem = matprop.get_Sigma_r(mat);
615 if (geom_type == HERMES_PLANAR)
617 result = D_elem[g] * int_grad_u_ext_grad_v<Real, Scalar>(n, wt, u_ext[g], v) +
618 Sigma_r_elem[g] * int_u_ext_v<Real, Scalar>(n, wt, u_ext[g], v);
622 if (geom_type == HERMES_AXISYM_X)
624 result = D_elem[g] * int_y_grad_u_ext_grad_v<Real, Scalar>(n, wt, u_ext[g], v, e) +
625 Sigma_r_elem[g] * int_y_u_ext_v<Real, Scalar>(n, wt, u_ext[g], v, e);
629 result = D_elem[g] * int_x_grad_u_ext_grad_v<Real, Scalar>(n, wt, u_ext[g], v, e) +
630 Sigma_r_elem[g] * int_x_u_ext_v<Real, Scalar>(n, wt, u_ext[g], v, e);
636 template<
typename ScalarClass>
637 template<
typename Real,
typename Scalar>
638 Scalar FissionYield::Jacobian<ScalarClass>::matrix_form(
int n,
double *wt, Func<Scalar> *u_ext[], Func<Real> *u,
639 Func<Real> *v, GeomVol<Real> *e, Func<Scalar> **ext)
const
641 if (!matprop.get_fission_multigroup_structure()[gto])
644 Scalar result = Scalar(0);
645 if (geom_type == HERMES_PLANAR) result = int_u_v<Real, Scalar>(n, wt, u, v);
648 if (geom_type == HERMES_AXISYM_X) result = int_y_u_v<Real, Scalar>(n, wt, u, v, e);
649 else result = int_x_u_v<Real, Scalar>(n, wt, u, v, e);
652 std::string mat = this->mesh->get_element_markers_conversion().get_user_marker(e->elem_marker).marker;
653 rank1 nu_elem = matprop.get_nu(mat);
654 rank1 Sigma_f_elem = matprop.get_Sigma_f(mat);
655 rank1 chi_elem = matprop.get_chi(mat);
657 return result * chi_elem[gto] * nu_elem[gfrom] * Sigma_f_elem[gfrom];
660 template<
typename ScalarClass>
661 template<
typename Real,
typename Scalar>
662 Scalar FissionYield::OuterIterationForm<ScalarClass>::vector_form(
int n,
double *wt, Func<Scalar> *u_ext[],
663 Func<Real> *v, GeomVol<Real> *e, Func<Scalar> **ext)
const
665 if (!matprop.get_fission_multigroup_structure()[g])
668 std::string mat = this->mesh->get_element_markers_conversion().get_user_marker(e->elem_marker).marker;
669 rank1 nu_elem = matprop.get_nu(mat);
670 rank1 Sigma_f_elem = matprop.get_Sigma_f(mat);
671 rank1 chi_elem = matprop.get_chi(mat);
673 Scalar result = Scalar(0);
674 for (
int i = 0; i < n; i++)
676 Scalar local_res = Scalar(0);
677 for (
int gfrom = 0; gfrom < this->wf->get_ext().size(); gfrom++)
678 local_res += nu_elem[gfrom] * Sigma_f_elem[gfrom] * ext[gfrom]->val[i];
680 local_res = local_res * wt[i] * v->val[i];
682 if (geom_type == HERMES_AXISYM_X)
683 local_res = local_res * e->y[i];
684 else if (geom_type == HERMES_AXISYM_Y)
685 local_res = local_res * e->x[i];
690 return result * chi_elem[g] / keff;
693 template<
typename ScalarClass>
694 template<
typename Real,
typename Scalar>
695 Scalar FissionYield::Residual<ScalarClass>::vector_form(
int n,
double *wt, Func<Scalar> *u_ext[],
696 Func<Real> *v, GeomVol<Real> *e, Func<Scalar> **ext)
const
698 if (!matprop.get_fission_multigroup_structure()[gto])
701 Scalar result = Scalar(0);
702 if (geom_type == HERMES_PLANAR) result = int_u_ext_v<Real, Scalar>(n, wt, u_ext[gfrom], v);
705 if (geom_type == HERMES_AXISYM_X) result = int_y_u_ext_v<Real, Scalar>(n, wt, u_ext[gfrom], v, e);
706 else result = int_x_u_ext_v<Real, Scalar>(n, wt, u_ext[gfrom], v, e);
709 std::string mat = this->mesh->get_element_markers_conversion().get_user_marker(e->elem_marker).marker;
710 rank1 nu_elem = matprop.get_nu(mat);
711 rank1 Sigma_f_elem = matprop.get_Sigma_f(mat);
712 rank1 chi_elem = matprop.get_chi(mat);
714 return result * chi_elem[gto] * nu_elem[gfrom] * Sigma_f_elem[gfrom];
717 template<
typename ScalarClass>
718 template<
typename Real,
typename Scalar>
719 Scalar Scattering::Jacobian<ScalarClass>::matrix_form(
int n,
double *wt, Func<Scalar> *u_ext[], Func<Real> *u,
720 Func<Real> *v, GeomVol<Real> *e, Func<Scalar> **ext)
const
722 Scalar result = Scalar(0);
723 if (geom_type == HERMES_PLANAR) result = int_u_v<Real, Scalar>(n, wt, u, v);
726 if (geom_type == HERMES_AXISYM_X) result = int_y_u_v<Real, Scalar>(n, wt, u, v, e);
727 else result = int_x_u_v<Real, Scalar>(n, wt, u, v, e);
730 return result * matprop.get_Sigma_s(this->mesh->get_element_markers_conversion().get_user_marker(e->elem_marker).marker)[gto][gfrom];
733 template<
typename ScalarClass>
734 template<
typename Real,
typename Scalar>
735 Scalar Scattering::Residual<ScalarClass>::vector_form(
int n,
double *wt, Func<Scalar> *u_ext[],
736 Func<Real> *v, GeomVol<Real> *e, Func<Scalar> **ext)
const
738 Scalar result = Scalar(0);
739 if (geom_type == HERMES_PLANAR) result = int_u_ext_v<Real, Scalar>(n, wt, u_ext[gfrom], v);
742 if (geom_type == HERMES_AXISYM_X) result = int_y_u_ext_v<Real, Scalar>(n, wt, u_ext[gfrom], v, e);
743 else result = int_x_u_ext_v<Real, Scalar>(n, wt, u_ext[gfrom], v, e);
746 return result * matprop.get_Sigma_s(this->mesh->get_element_markers_conversion().get_user_marker(e->elem_marker).marker)[gto][gfrom];
749 template<
typename ScalarClass>
750 template<
typename Real,
typename Scalar>
751 Scalar ExternalSources::LinearForm<ScalarClass>::vector_form(
int n,
double *wt, Func<Scalar> *u_ext[],
752 Func<Real> *v, GeomVol<Real> *e, Func<Scalar> **ext)
const
754 std::string mat = this->mesh->get_element_markers_conversion().get_user_marker(e->elem_marker).marker;
756 if (geom_type == HERMES_PLANAR)
757 return matprop.get_src(mat)[g] * int_v<Real>(n, wt, v);
760 if (geom_type == HERMES_AXISYM_X)
761 return matprop.get_src(mat)[g] * int_y_v<Real>(n, wt, v, e);
763 return matprop.get_src(mat)[g] * int_x_v<Real>(n, wt, v, e);
769 namespace CompleteWeakForms
773 template<
typename Scalar>
774 void DefaultWeakFormFixedSource<Scalar>::lhs_init(
unsigned int G,
const MaterialPropertyMaps& matprop, MeshSharedPtr mesh,
777 bool2 Ss_nnz = matprop.get_scattering_multigroup_structure();
778 bool1 chi_nnz = matprop.get_fission_multigroup_structure();
780 for (
unsigned int gto = 0; gto < G; gto++)
782 this->add_matrix_form(
new DiffusionReaction::Jacobian<Scalar>(gto, matprop, mesh, geom_type));
783 this->add_vector_form(
new DiffusionReaction::Residual<Scalar>(gto, matprop, mesh, geom_type));
785 for (
unsigned int gfrom = 0; gfrom < G; gfrom++)
787 if (Ss_nnz[gto][gfrom])
789 this->add_matrix_form(
new Scattering::Jacobian<Scalar>(gto, gfrom, matprop, mesh, geom_type));
790 this->add_vector_form(
new Scattering::Residual<Scalar>(gto, gfrom, matprop, mesh, geom_type));
795 this->add_matrix_form(
new FissionYield::Jacobian<Scalar>(gto, gfrom, matprop, mesh, geom_type));
796 this->add_vector_form(
new FissionYield::Residual<Scalar>(gto, gfrom, matprop, mesh, geom_type));
802 template<
typename Scalar>
803 DefaultWeakFormFixedSource<Scalar>::DefaultWeakFormFixedSource(
const MaterialPropertyMaps& matprop, MeshSharedPtr mesh,
804 GeomType geom_type) : WeakForm<Scalar>(matprop.get_G())
806 lhs_init(matprop.get_G(), matprop, mesh, geom_type);
807 for (
unsigned int gto = 0; gto < matprop.get_G(); gto++)
808 this->add_vector_form(
new ExternalSources::LinearForm<Scalar>(gto, matprop, mesh, geom_type));
811 template<
typename Scalar>
812 DefaultWeakFormFixedSource<Scalar>::DefaultWeakFormFixedSource(
const MaterialPropertyMaps& matprop, MeshSharedPtr mesh,
813 Hermes2DFunction<Scalar>*f_src,
std::string src_area,
814 GeomType geom_type) : WeakForm<Scalar>(matprop.get_G())
816 lhs_init(matprop.get_G(), matprop, mesh, geom_type);
817 for (
unsigned int gto = 0; gto < matprop.get_G(); gto++)
818 this->add_vector_form(
new WeakFormsH1::DefaultVectorFormVol<Scalar>(gto, src_area, f_src, geom_type));
821 template<
typename Scalar>
822 DefaultWeakFormFixedSource<Scalar>::DefaultWeakFormFixedSource(
const MaterialPropertyMaps& matprop, MeshSharedPtr mesh,
823 Hermes2DFunction<Scalar>*f_src,
824 std::vector<std::string> src_areas,
825 GeomType geom_type) : WeakForm<Scalar>(matprop.get_G())
827 lhs_init(matprop.get_G(), matprop, mesh, geom_type);
828 for (
unsigned int gto = 0; gto < matprop.get_G(); gto++)
829 this->add_vector_form(
new WeakFormsH1::DefaultVectorFormVol<Scalar>(gto, src_areas, f_src, geom_type));
832 template<
typename Scalar>
833 DefaultWeakFormFixedSource<Scalar>::DefaultWeakFormFixedSource(
const MaterialPropertyMaps& matprop, MeshSharedPtr mesh,
834 const std::vector<Hermes2DFunction<Scalar>*>& f_src,
836 GeomType geom_type) : WeakForm<Scalar>(matprop.get_G())
838 if (f_src.size() != matprop.get_G())
839 throw Hermes::Exceptions::Exception(E_INVALID_SIZE);
841 lhs_init(matprop.get_G(), matprop, mesh, geom_type);
842 for (
unsigned int gto = 0; gto < matprop.get_G(); gto++)
843 this->add_vector_form(
new WeakFormsH1::DefaultVectorFormVol<Scalar>(gto, src_area, f_src[gto], geom_type));
846 template<
typename Scalar>
847 DefaultWeakFormFixedSource<Scalar>::DefaultWeakFormFixedSource(
const MaterialPropertyMaps& matprop, MeshSharedPtr mesh,
848 const std::vector<Hermes2DFunction<Scalar>*>& f_src,
849 std::vector<std::string> src_areas,
850 GeomType geom_type) : WeakForm<Scalar>(matprop.get_G())
852 if (f_src.size() != matprop.get_G())
853 throw Hermes::Exceptions::Exception(E_INVALID_SIZE);
855 lhs_init(matprop.get_G(), matprop, mesh, geom_type);
856 for (
unsigned int gto = 0; gto < matprop.get_G(); gto++)
857 this->add_vector_form(
new WeakFormsH1::DefaultVectorFormVol<Scalar>(gto, src_areas, f_src[gto], geom_type));
860 template<
typename Scalar>
861 DefaultWeakFormSourceIteration<Scalar>::DefaultWeakFormSourceIteration(
const MaterialPropertyMaps& matprop, MeshSharedPtr mesh,
862 std::vector<MeshFunctionSharedPtr<Scalar> >& iterates,
863 double initial_keff_guess,
864 GeomType geom_type) : WeakForm<Scalar>(matprop.get_G())
866 bool2 Ss_nnz = matprop.get_scattering_multigroup_structure();
868 for (
unsigned int gto = 0; gto < matprop.get_G(); gto++)
870 this->add_matrix_form(
new DiffusionReaction::Jacobian<Scalar>(gto, matprop, mesh, geom_type));
871 this->add_vector_form(
new DiffusionReaction::Residual<Scalar>(gto, matprop, mesh, geom_type));
873 for (
unsigned int gfrom = 0; gfrom < matprop.get_G(); gfrom++)
875 if (Ss_nnz[gto][gfrom])
877 this->add_matrix_form(
new Scattering::Jacobian<Scalar>(gto, gfrom, matprop, mesh, geom_type));
878 this->add_vector_form(
new Scattering::Residual<Scalar>(gto, gfrom, matprop, mesh, geom_type));
882 FissionYield::OuterIterationForm<Scalar>* keff_iteration_form =
new FissionYield::OuterIterationForm<Scalar>(gto, matprop, mesh, iterates, initial_keff_guess, geom_type);
884 keff_iteration_forms.push_back(keff_iteration_form);
885 this->add_vector_form(keff_iteration_form);
889 template<
typename Scalar>
890 void DefaultWeakFormSourceIteration<Scalar>::update_keff(
double new_keff)
897 for (
int i = 0; i < keff_iteration_forms.size(); i++)
898 keff_iteration_forms[i]->update_keff(new_keff);
903 namespace SupportClasses
905 void SourceFilter::filter_fn(
int n,
const std::vector<const double*>& values,
double* result)
907 for (
int i = 0; i < n; i++)
910 for (
unsigned int j = 0; j < values.size(); j++)
911 result[i] += nu[j] * Sigma_f[j] * values.at(j)[i];
916 namespace Monoenergetic
920 template class HERMES_API DefaultWeakFormFixedSource < double > ;
921 template class HERMES_API DefaultWeakFormFixedSource < std::complex<double> > ;
926 namespace ElementaryForms
930 template class HERMES_API FissionYield::Jacobian < double > ;
931 template class HERMES_API FissionYield::Jacobian < std::complex<double> > ;
932 template class HERMES_API FissionYield::OuterIterationForm < double > ;
933 template class HERMES_API FissionYield::OuterIterationForm < std::complex<double> > ;
934 template class HERMES_API FissionYield::Residual < double > ;
935 template class HERMES_API FissionYield::Residual < std::complex<double> > ;
937 template class HERMES_API VacuumBoundaryCondition::Jacobian < double > ;
938 template class HERMES_API VacuumBoundaryCondition::Jacobian < std::complex<double> > ;
939 template class HERMES_API VacuumBoundaryCondition::Residual < double > ;
940 template class HERMES_API VacuumBoundaryCondition::Residual < std::complex<double> > ;
942 template double VacuumBoundaryCondition::Jacobian<double>::matrix_form<double,
double>(
int n,
double *wt, Func<double> *u_ext[], Func<double> *u,
943 Func<double> *v, GeomSurf<double> *e, Func<double> **ext)
const;
945 template double VacuumBoundaryCondition::Residual<double>::vector_form<double,
double>(
int n,
double *wt, Func<double> *u_ext[],
946 Func<double> *v, GeomSurf<double> *e, Func<double> **ext)
const;
948 template double DiffusionReaction::Jacobian<double>::matrix_form<double,
double>(
int n,
double *wt, Func<double> *u_ext[], Func<double> *u,
949 Func<double> *v, GeomVol<double> *e, Func<double> **ext)
const;
951 template double DiffusionReaction::Residual<double>::vector_form<double,
double>(
int n,
double *wt, Func<double> *u_ext[],
952 Func<double> *v, GeomVol<double> *e, Func<double> **ext)
const;
954 template double FissionYield::Jacobian<double>::matrix_form<double,
double>(
int n,
double *wt, Func<double> *u_ext[], Func<double> *u,
955 Func<double> *v, GeomVol<double> *e, Func<double> **ext)
const;
957 template double FissionYield::OuterIterationForm<double>::vector_form<double,
double>(
int n,
double *wt, Func<double> *u_ext[],
958 Func<double> *v, GeomVol<double> *e, Func<double> **ext)
const;
960 template double FissionYield::Residual<double>::vector_form<double,
double>(
int n,
double *wt, Func<double> *u_ext[],
961 Func<double> *v, GeomVol<double> *e, Func<double> **ext)
const;
963 template double Scattering::Jacobian<double>::matrix_form<double,
double>(
int n,
double *wt, Func<double> *u_ext[], Func<double> *u,
964 Func<double> *v, GeomVol<double> *e, Func<double> **ext)
const;
966 template double Scattering::Residual<double>::vector_form<double,
double>(
int n,
double *wt, Func<double> *u_ext[],
967 Func<double> *v, GeomVol<double> *e, Func<double> **ext)
const;
969 template double ExternalSources::LinearForm<double>::vector_form<double,
double>(
int n,
double *wt, Func<double> *u_ext[],
970 Func<double> *v, GeomVol<double> *e, Func<double> **ext)
const;
974 namespace CompleteWeakForms
978 template class HERMES_API DefaultWeakFormFixedSource < double > ;
979 template class HERMES_API DefaultWeakFormFixedSource < std::complex<double> > ;
981 template class HERMES_API DefaultWeakFormSourceIteration < double > ;
982 template class HERMES_API DefaultWeakFormSourceIteration < std::complex<double> > ;
::xsd::cxx::tree::string< char, simple_type > string
C++ type corresponding to the string XML Schema built-in type.
GeomType
Geometrical type of weak forms.