16 #include "weakforms_neutronics.h"
24 namespace WeakFormsNeutronics
26 namespace Monoenergetic
30 template<
typename Scalar>
31 DefaultWeakFormFixedSource<Scalar>::DefaultWeakFormFixedSource( Hermes::vector<std::string> regions,
32 Hermes::vector<double> D_map,
33 Hermes::vector<double> Sigma_a_map,
34 Hermes::vector<double> Q_map ) : WeakForm<Scalar>(1)
36 using namespace WeakFormsH1;
38 for (
unsigned int i = 0; i < regions.size(); i++)
42 this->add_matrix_form(
new DefaultJacobianDiffusion<Scalar>(0, 0, regions[i],
new Hermes1DFunction<Scalar>(D_map[i]),
45 this->add_matrix_form(
new DefaultMatrixFormVol<Scalar>(0, 0, regions[i],
new Hermes2DFunction<Scalar>(Sigma_a_map[i]),
50 this->add_vector_form(
new DefaultResidualDiffusion<Scalar>(0, regions[i],
new Hermes1DFunction<Scalar>(D_map[i])));
52 this->add_vector_form(
new DefaultResidualVol<Scalar>(0, regions[i],
new Hermes2DFunction<Scalar>(Sigma_a_map[i])));
54 this->add_vector_form(
new DefaultVectorFormVol<Scalar>(0, regions[i],
new Hermes2DFunction<Scalar>(-Q_map[i])));
62 namespace MaterialProperties
66 void MaterialPropertyMaps::extend_to_multigroup(
const MaterialPropertyMap0& mrsg_map,
67 MaterialPropertyMap1 *mrmg_map)
70 this->warn(W_MG_EXTENSION);
72 MaterialPropertyMap0::const_iterator it;
73 for (it = mrsg_map.begin(); it != mrsg_map.end(); ++it)
74 (*mrmg_map)[it->first].assign(G, it->second);
77 void MaterialPropertyMaps::extend_to_multiregion(
const rank1& srmg_array,
78 MaterialPropertyMap1 *mrmg_map)
80 if(materials_list.empty())
81 throw Hermes::Exceptions::Exception(E_MR_EXTENSION);
83 std::set<std::string>::const_iterator it;
84 for (it = materials_list.begin(); it != materials_list.end(); ++it)
85 (*mrmg_map)[*it] = srmg_array;
88 void MaterialPropertyMaps::extend_to_multiregion_multigroup(
const rank0& srsg_value,
89 MaterialPropertyMap1 *mrmg_map)
91 if(materials_list.empty())
92 throw Hermes::Exceptions::Exception(E_MR_EXTENSION);
94 this->warn(W_MG_EXTENSION);
96 std::set<std::string>::const_iterator it;
97 for (it = materials_list.begin(); it != materials_list.end(); ++it)
98 (*mrmg_map)[*it].assign(G, srsg_value);
101 void MaterialPropertyMaps::fill_with(
double c, MaterialPropertyMap1 *mrmg_map)
103 if(materials_list.empty())
104 throw Hermes::Exceptions::Exception(E_MR_EXTENSION);
106 std::set<std::string>::const_iterator it;
107 for (it = materials_list.begin(); it != materials_list.end(); ++it)
108 (*mrmg_map)[*it].assign(G, c);
111 void MaterialPropertyMaps::validate()
113 using namespace ValidationFunctors;
115 if(fission_multigroup_structure.empty())
116 fission_multigroup_structure = bool1(G,
true);
120 fill_with(0.0, &chi);
121 MaterialPropertyMap1::iterator it = chi.begin();
122 for ( ; it != chi.end(); ++it)
124 fission_multigroup_structure = bool1(G,
false);
125 fission_multigroup_structure[0] =
true;
128 if(nu.empty() && !nuSigma_f.empty() && !Sigma_f.empty())
129 nu = NDArrayMapOp::divide<rank1>(nuSigma_f, Sigma_f);
130 else if(nuSigma_f.empty() && !nu.empty() && !Sigma_f.empty())
131 nuSigma_f = NDArrayMapOp::multiply<rank1>(nu, Sigma_f);
132 else if(Sigma_f.empty() && !nuSigma_f.empty() && !nu.empty())
133 Sigma_f = NDArrayMapOp::divide<rank1>(nuSigma_f, nu);
134 else if(!Sigma_f.empty() && !nuSigma_f.empty() && !nu.empty())
136 MaterialPropertyMap1 diff = NDArrayMapOp::subtract<rank1>(nuSigma_f,
137 NDArrayMapOp::multiply<rank1>(nu, Sigma_f) );
138 std::for_each(diff.begin(), diff.end(), ensure_trivial());
142 this->warn(W_NO_FISSION);
144 fill_with(0.0, &chi);
145 fill_with(0.0, &Sigma_f);
148 if((nu.size() != Sigma_f.size()) || (nu.size() != chi.size()))
149 throw Hermes::Exceptions::Exception(E_NONMATCHING_PROPERTIES);
151 if(Sigma_f.size() > 0)
153 std::for_each(nu.begin(), nu.end(), ensure_size(G));
154 std::for_each(Sigma_f.begin(), Sigma_f.end(), ensure_size(G));
155 std::for_each(chi.begin(), chi.end(), ensure_size(G));
158 if(Sigma_a.size() > 0)
163 MaterialPropertyMap1::const_iterator ita = Sigma_a.begin();
164 MaterialPropertyMap1::const_iterator itf = Sigma_f.begin();
165 for ( ; ita != Sigma_a.end(); ++ita, ++itf)
167 rank1::const_iterator a = ita->second.begin();
168 rank1::const_iterator f = itf->second.begin();
170 for ( ; a != ita->second.end(); ++a, ++f)
172 this->warn(W_SA_LT_SF);
177 const rank1& MaterialPropertyMaps::get_Sigma_f(
std::string material)
const
179 if(material ==
"-999")
return this->Sigma_f.begin()->second;
183 MaterialPropertyMap1::const_iterator data = this->Sigma_f.find(material);
184 if(data != this->Sigma_f.end())
188 throw Hermes::Exceptions::Exception(E_INVALID_MARKER);
189 return *(
new rank1());
192 const rank1& MaterialPropertyMaps::get_nu(
std::string material)
const
194 if(material ==
"-999")
return this->nu.begin()->second;
196 MaterialPropertyMap1::const_iterator data = this->nu.find(material);
197 if(data != this->nu.end())
201 throw Hermes::Exceptions::Exception(E_INVALID_MARKER);
202 return *(
new rank1());
205 const rank1& MaterialPropertyMaps::get_chi(
std::string material)
const
207 if(material ==
"-999")
return this->chi.begin()->second;
209 MaterialPropertyMap1::const_iterator data = this->chi.find(material);
210 if(data != this->chi.end())
214 throw Hermes::Exceptions::Exception(E_INVALID_MARKER);
215 return *(
new rank1());
219 std::ostream & operator<< (std::ostream& os,
const MaterialPropertyMaps& matprop)
224 os << setw(12) <<
"target group" << setw(10) <<
"chi" << setw(10) <<
"nu";
225 os << setw(10) <<
"Sigma_f" << std::endl;
227 MaterialPropertyMap1::const_iterator data_elem = matprop.chi.begin();
228 for ( ; data_elem != matprop.chi.end(); ++data_elem)
230 string mat = data_elem->first;
232 os << setw(80) << setfill(
'-') <<
' ' << std::endl << setfill(
' ');
233 os << setw(40) << mat << std::endl;
234 os << setw(80) << setfill(
'-') <<
' ' << std::endl << setfill(
' ');
235 for (
unsigned int gto = 0; gto < matprop.G; gto++)
237 os << setw(6) << gto << setw(6) <<
' ';
238 os << setw(10) << matprop.get_chi(mat)[gto];
239 os << setw(10) << matprop.get_nu(mat)[gto];
240 os << setw(10) << matprop.get_Sigma_f(mat)[gto];
253 MaterialPropertyMap1 MaterialPropertyMaps::extract_map2_diagonals(
const MaterialPropertyMap2& map2)
255 MaterialPropertyMap1 diags;
257 MaterialPropertyMap2::const_iterator map2_it = map2.begin();
258 for ( ; map2_it != map2.end(); ++map2_it)
260 diags[map2_it->first].reserve(G);
261 for (
unsigned int g = 0; g < G; g++)
262 diags[map2_it->first].push_back(map2_it->second[g][g]);
268 MaterialPropertyMap1 MaterialPropertyMaps::sum_map2_columns(
const MaterialPropertyMap2& map2)
270 MaterialPropertyMap1 summed;
272 MaterialPropertyMap2::const_iterator map2_it = map2.begin();
273 for ( ; map2_it != map2.end(); ++map2_it)
275 summed[map2_it->first].reserve(G);
276 for (
unsigned int gfrom = 0; gfrom < G; gfrom++)
280 for (
unsigned int gto = 0; gto < G; gto++)
281 sum += map2_it->second[gto][gfrom];
283 summed[map2_it->first].push_back(sum);
290 MaterialPropertyMap2 MaterialPropertyMaps::create_map2_by_diagonals(
const MaterialPropertyMap1& diags)
292 MaterialPropertyMap2 map2;
294 MaterialPropertyMap1::const_iterator diags_it = diags.begin();
295 for ( ; diags_it != diags.end(); ++diags_it)
297 map2[diags_it->first].resize(G, rank1(G, 0.0));
299 for (
unsigned int g = 0; g < G; g++)
300 map2[diags_it->first][g][g] = diags_it->second[g];
306 void MaterialPropertyMaps::fill_with(
double c, MaterialPropertyMap2 *mrmg_map)
308 std::set<std::string>::const_iterator it;
309 for (it = materials_list.begin(); it != materials_list.end(); ++it)
310 (*mrmg_map)[*it].assign(G, rank1(G, c));
313 void MaterialPropertyMaps::validate()
315 Common::MaterialPropertyMaps::validate();
317 bool D_given = !D.empty();
318 bool Sigma_r_given = !Sigma_r.empty();
319 bool Sigma_s_given = !Sigma_s.empty();
320 bool Sigma_t_given = !Sigma_t.empty();
321 bool Sigma_a_given = !Sigma_a.empty();
322 bool Sigma_f_given = !Sigma_f.empty();
323 bool src_given = !src.empty();
336 Sigma_s = create_map2_by_diagonals(Common::NDArrayMapOp::subtract<rank1>(Sigma_t, Sigma_a));
342 this->warn(W_NO_SCATTERING);
343 fill_with(0.0, &Sigma_s);
346 Sigma_s_given =
true;
355 this->warn(W_NO_SCATTERING);
356 fill_with(0.0, &Sigma_s);
357 Sigma_s_given =
true;
361 Sigma_t = Common::NDArrayMapOp::add<rank1>(Sigma_a, sum_map2_columns(Sigma_s));
368 throw Hermes::Exceptions::Exception(E_INSUFFICIENT_DATA);
371 Sigma_t_given =
true;
374 Sigma_r = Common::NDArrayMapOp::subtract<rank1>(Sigma_t, extract_map2_diagonals(Sigma_s));
375 Sigma_r_given =
true;
380 if(scattering_multigroup_structure.empty())
381 scattering_multigroup_structure = bool2(G, Hermes::vector<bool>(G,
true));
391 Sigma_s = create_map2_by_diagonals(Common::NDArrayMapOp::subtract<rank1>(Sigma_t, Sigma_r));
393 scattering_multigroup_structure = bool2(G, Hermes::vector<bool>(G,
false));
394 for (
unsigned int gto = 0; gto < G; gto++)
395 for (
unsigned int gfrom = 0; gfrom < G; gfrom++)
397 scattering_multigroup_structure[gto][gfrom] =
true;
401 this->warn(W_NO_SCATTERING);
402 fill_with(0.0, &Sigma_s);
403 scattering_multigroup_structure = bool2(G, Hermes::vector<bool>(G,
false));
406 Sigma_s_given =
true;
413 MaterialPropertyMap1::const_iterator Sr_elem = Sigma_r.begin();
414 for ( ; Sr_elem != Sigma_r.end(); ++Sr_elem)
415 for (
unsigned int g = 0; g < G; g++)
416 D[Sr_elem->first][g] = 1./(3.*Sr_elem->second[g]);
421 if((D.size() != Sigma_r.size()) || (D.size() != Sigma_s.size()) || (src_given && D.size() != src.size()))
422 throw Hermes::Exceptions::Exception(E_NONMATCHING_PROPERTIES);
424 using ValidationFunctors::ensure_size;
425 std::for_each(Sigma_s.begin(), Sigma_s.end(), ensure_size(G, G));
426 std::for_each(Sigma_r.begin(), Sigma_r.end(), ensure_size(G));
427 std::for_each(src.begin(), src.end(), ensure_size(G));
428 std::for_each(D.begin(), D.end(), ensure_size(G));
431 const rank2& MaterialPropertyMaps::get_Sigma_s(
std::string material)
const
433 if(material ==
"-999")
return this->Sigma_s.begin()->second;
437 MaterialPropertyMap2::const_iterator data = this->Sigma_s.find(material);
438 if(data != this->Sigma_s.end())
442 throw Hermes::Exceptions::Exception(E_INVALID_MARKER);
443 return *(
new rank2());
446 const rank1& MaterialPropertyMaps::get_Sigma_r(
std::string material)
const
448 if(material ==
"-999")
return this->Sigma_r.begin()->second;
450 MaterialPropertyMap1::const_iterator data = this->Sigma_r.find(material);
451 if(data != this->Sigma_r.end())
455 throw Hermes::Exceptions::Exception(E_INVALID_MARKER);
456 return *(
new rank1());
459 const rank1& MaterialPropertyMaps::get_D(
std::string material)
const
461 if(material ==
"-999")
return this->D.begin()->second;
463 MaterialPropertyMap1::const_iterator data = this->D.find(material);
464 if(data != this->D.end())
468 throw Hermes::Exceptions::Exception(E_INVALID_MARKER);
469 return *(
new rank1());
472 const rank1& MaterialPropertyMaps::get_src(
std::string material)
const
474 if(material ==
"-999")
return this->src.begin()->second;
476 MaterialPropertyMap1::const_iterator data = this->src.find(material);
477 if(data != this->src.end())
481 throw Hermes::Exceptions::Exception(E_INVALID_MARKER);
482 return *(
new rank1());
486 HERMES_API std::ostream & operator<< (std::ostream& os,
const MaterialPropertyMaps& matprop)
490 os << static_cast<const Common::MaterialPropertyMaps&>(matprop) << std::endl;
492 os << setw(12) <<
"target group" << setw(10) <<
"D" << setw(10) <<
"Sigma_r";
493 os << setw(10) <<
"ext. src" << setw(22) <<
"Sigma_s" << std::endl;
495 MaterialPropertyMap1::const_iterator data_elem = matprop.Sigma_r.begin();
496 for ( ; data_elem != matprop.Sigma_r.end(); ++data_elem)
498 string mat = data_elem->first;
500 os << setw(80) << setfill(
'-') <<
' ' << std::endl << setfill(
' ');
501 os << setw(40) << mat << std::endl;
502 os << setw(80) << setfill(
'-') <<
' ' << std::endl << setfill(
' ');
503 for (
unsigned int gto = 0; gto < matprop.G; gto++)
505 os << setw(6) << gto << setw(6) <<
' ';
506 os << setw(10) << matprop.get_D(mat)[gto];
507 os << setw(10) << matprop.get_Sigma_r(mat)[gto];
509 if(matprop.src.empty())
512 os << matprop.get_src(mat)[gto];
514 for (
unsigned int gfrom = 0; gfrom < matprop.G; gfrom++)
515 os << setw(8) << matprop.get_Sigma_s(mat)[gto][gfrom];
521 return os << std::endl;
526 namespace ElementaryForms
530 template<
typename ScalarClass>
531 template<
typename Real,
typename Scalar>
532 Scalar VacuumBoundaryCondition::Jacobian<ScalarClass>::matrix_form(
int n,
double *wt, Func<Scalar> *u_ext[], Func<Real> *u,
533 Func<Real> *v, Geom<Real> *e, Func<Scalar> **ext)
const
537 if(geom_type == HERMES_PLANAR)
538 result = 0.5 * int_u_v<Real, Scalar>(n, wt, u, v);
539 else if(geom_type == HERMES_AXISYM_X)
540 result = 0.5 * int_y_u_v<Real, Scalar>(n, wt, u, v, e);
542 result = 0.5 * int_x_u_v<Real, Scalar>(n, wt, u, v, e);
547 template<
typename ScalarClass>
548 template<
typename Real,
typename Scalar>
549 Scalar VacuumBoundaryCondition::Residual<ScalarClass>::vector_form(
int n,
double *wt, Func<Scalar> *u_ext[],
550 Func<Real> *v, Geom<Real> *e, Func<Scalar> **ext)
const
554 if(geom_type == HERMES_PLANAR)
555 result = 0.5 * int_u_ext_v<Real, Scalar>(n, wt, u_ext[g], v);
556 else if(geom_type == HERMES_AXISYM_X)
557 result = 0.5 * int_y_u_ext_v<Real, Scalar>(n, wt, u_ext[g], v, e);
559 result = 0.5 * int_x_u_ext_v<Real, Scalar>(n, wt, u_ext[g], v, e);
564 template<
typename ScalarClass>
565 template<
typename Real,
typename Scalar>
566 Scalar DiffusionReaction::Jacobian<ScalarClass>::matrix_form(
int n,
double *wt, Func<Scalar> *u_ext[], Func<Real> *u,
567 Func<Real> *v, Geom<Real> *e, Func<Scalar> **ext)
const
571 std::string mat = this->mesh->get_element_markers_conversion().get_user_marker(e->elem_marker).marker;
572 rank1 D_elem = matprop.get_D(mat);
573 rank1 Sigma_r_elem = matprop.get_Sigma_r(mat);
575 if(geom_type == HERMES_PLANAR)
577 result = D_elem[g] * int_grad_u_grad_v<Real, Scalar>(n, wt, u, v) +
578 Sigma_r_elem[g] * int_u_v<Real, Scalar>(n, wt, u, v);
582 if(geom_type == HERMES_AXISYM_X)
584 result = D_elem[g] * int_y_grad_u_grad_v<Real, Scalar>(n, wt, u, v, e) +
585 Sigma_r_elem[g] * int_y_u_v<Real, Scalar>(n, wt, u, v, e);
589 result = D_elem[g] * int_x_grad_u_grad_v<Real, Scalar>(n, wt, u, v, e) +
590 Sigma_r_elem[g] * int_x_u_v<Real, Scalar>(n, wt, u, v, e);
596 template<
typename ScalarClass>
597 template<
typename Real,
typename Scalar>
598 Scalar DiffusionReaction::Residual<ScalarClass>::vector_form(
int n,
double *wt, Func<Scalar> *u_ext[],
599 Func<Real> *v, Geom<Real> *e, Func<Scalar> **ext)
const
603 std::string mat = this->mesh->get_element_markers_conversion().get_user_marker(e->elem_marker).marker;
604 rank1 D_elem = matprop.get_D(mat);
605 rank1 Sigma_r_elem = matprop.get_Sigma_r(mat);
607 if(geom_type == HERMES_PLANAR)
609 result = D_elem[g] * int_grad_u_ext_grad_v<Real, Scalar>(n, wt, u_ext[g], v) +
610 Sigma_r_elem[g] * int_u_ext_v<Real, Scalar>(n, wt, u_ext[g], v);
614 if(geom_type == HERMES_AXISYM_X)
616 result = D_elem[g] * int_y_grad_u_ext_grad_v<Real, Scalar>(n, wt, u_ext[g], v, e) +
617 Sigma_r_elem[g] * int_y_u_ext_v<Real, Scalar>(n, wt, u_ext[g], v, e);
621 result = D_elem[g] * int_x_grad_u_ext_grad_v<Real, Scalar>(n, wt, u_ext[g], v, e) +
622 Sigma_r_elem[g] * int_x_u_ext_v<Real, Scalar>(n, wt, u_ext[g], v, e);
628 template<
typename ScalarClass>
629 template<
typename Real,
typename Scalar>
630 Scalar FissionYield::Jacobian<ScalarClass>::matrix_form(
int n,
double *wt, Func<Scalar> *u_ext[], Func<Real> *u,
631 Func<Real> *v, Geom<Real> *e, Func<Scalar> **ext)
const
633 if(!matprop.get_fission_multigroup_structure()[gto])
636 Scalar result = Scalar(0);
637 if(geom_type == HERMES_PLANAR) result = int_u_v<Real, Scalar>(n, wt, u, v);
640 if(geom_type == HERMES_AXISYM_X) result = int_y_u_v<Real, Scalar>(n, wt, u, v, e);
641 else result = int_x_u_v<Real, Scalar>(n, wt, u, v, e);
644 std::string mat = this->mesh->get_element_markers_conversion().get_user_marker(e->elem_marker).marker;
645 rank1 nu_elem = matprop.get_nu(mat);
646 rank1 Sigma_f_elem = matprop.get_Sigma_f(mat);
647 rank1 chi_elem = matprop.get_chi(mat);
649 return result * chi_elem[gto] * nu_elem[gfrom] * Sigma_f_elem[gfrom];
652 template<
typename ScalarClass>
653 template<
typename Real,
typename Scalar>
654 Scalar FissionYield::OuterIterationForm<ScalarClass>::vector_form(
int n,
double *wt, Func<Scalar> *u_ext[],
655 Func<Real> *v, Geom<Real> *e, Func<Scalar> **ext)
const
657 if(!matprop.get_fission_multigroup_structure()[g])
660 std::string mat = this->mesh->get_element_markers_conversion().get_user_marker(e->elem_marker).marker;
661 rank1 nu_elem = matprop.get_nu(mat);
662 rank1 Sigma_f_elem = matprop.get_Sigma_f(mat);
663 rank1 chi_elem = matprop.get_chi(mat);
665 Scalar result = Scalar(0);
666 for (
int i = 0; i < n; i++)
668 Scalar local_res = Scalar(0);
669 for (
int gfrom = 0; gfrom < this->wf->get_ext().size(); gfrom++)
670 local_res += nu_elem[gfrom] * Sigma_f_elem[gfrom] * ext[gfrom]->val[i];
672 local_res = local_res * wt[i] * v->val[i];
674 if(geom_type == HERMES_AXISYM_X)
675 local_res = local_res * e->y[i];
676 else if(geom_type == HERMES_AXISYM_Y)
677 local_res = local_res * e->x[i];
682 return result * chi_elem[g] / keff;
685 template<
typename ScalarClass>
686 template<
typename Real,
typename Scalar>
687 Scalar FissionYield::Residual<ScalarClass>::vector_form(
int n,
double *wt, Func<Scalar> *u_ext[],
688 Func<Real> *v, Geom<Real> *e, Func<Scalar> **ext )
const
690 if(!matprop.get_fission_multigroup_structure()[gto])
693 Scalar result = Scalar(0);
694 if(geom_type == HERMES_PLANAR) result = int_u_ext_v<Real, Scalar>(n, wt, u_ext[gfrom], v);
697 if(geom_type == HERMES_AXISYM_X) result = int_y_u_ext_v<Real, Scalar>(n, wt, u_ext[gfrom], v, e);
698 else result = int_x_u_ext_v<Real, Scalar>(n, wt, u_ext[gfrom], v, e);
701 std::string mat = this->mesh->get_element_markers_conversion().get_user_marker(e->elem_marker).marker;
702 rank1 nu_elem = matprop.get_nu(mat);
703 rank1 Sigma_f_elem = matprop.get_Sigma_f(mat);
704 rank1 chi_elem = matprop.get_chi(mat);
706 return result * chi_elem[gto] * nu_elem[gfrom] * Sigma_f_elem[gfrom];
709 template<
typename ScalarClass>
710 template<
typename Real,
typename Scalar>
711 Scalar Scattering::Jacobian<ScalarClass>::matrix_form(
int n,
double *wt, Func<Scalar> *u_ext[], Func<Real> *u,
712 Func<Real> *v, Geom<Real> *e, Func<Scalar> **ext )
const
714 Scalar result = Scalar(0);
715 if(geom_type == HERMES_PLANAR) result = int_u_v<Real, Scalar>(n, wt, u, v);
718 if(geom_type == HERMES_AXISYM_X) result = int_y_u_v<Real, Scalar>(n, wt, u, v, e);
719 else result = int_x_u_v<Real, Scalar>(n, wt, u, v, e);
722 return result * matprop.get_Sigma_s(this->mesh->get_element_markers_conversion().get_user_marker(e->elem_marker).marker)[gto][gfrom];
725 template<
typename ScalarClass>
726 template<
typename Real,
typename Scalar>
727 Scalar Scattering::Residual<ScalarClass>::vector_form(
int n,
double *wt, Func<Scalar> *u_ext[],
728 Func<Real> *v, Geom<Real> *e, Func<Scalar> **ext )
const
730 Scalar result = Scalar(0);
731 if(geom_type == HERMES_PLANAR) result = int_u_ext_v<Real, Scalar>(n, wt, u_ext[gfrom], v);
734 if(geom_type == HERMES_AXISYM_X) result = int_y_u_ext_v<Real, Scalar>(n, wt, u_ext[gfrom], v, e);
735 else result = int_x_u_ext_v<Real, Scalar>(n, wt, u_ext[gfrom], v, e);
738 return result * matprop.get_Sigma_s(this->mesh->get_element_markers_conversion().get_user_marker(e->elem_marker).marker)[gto][gfrom];
741 template<
typename ScalarClass>
742 template<
typename Real,
typename Scalar>
743 Scalar ExternalSources::LinearForm<ScalarClass>::vector_form(
int n,
double *wt, Func<Scalar> *u_ext[],
744 Func<Real> *v, Geom<Real> *e, Func<Scalar> **ext)
const
746 std::string mat = this->mesh->get_element_markers_conversion().get_user_marker(e->elem_marker).marker;
748 if(geom_type == HERMES_PLANAR)
749 return matprop.get_src(mat)[g] * int_v<Real>(n, wt, v);
752 if(geom_type == HERMES_AXISYM_X)
753 return matprop.get_src(mat)[g] * int_y_v<Real>(n, wt, v, e);
755 return matprop.get_src(mat)[g] * int_x_v<Real>(n, wt, v, e);
761 namespace CompleteWeakForms
765 template<
typename Scalar>
766 void DefaultWeakFormFixedSource<Scalar>::lhs_init(
unsigned int G,
const MaterialPropertyMaps& matprop, Mesh *mesh,
769 bool2 Ss_nnz = matprop.get_scattering_multigroup_structure();
770 bool1 chi_nnz = matprop.get_fission_multigroup_structure();
772 for (
unsigned int gto = 0; gto < G; gto++)
774 this->add_matrix_form(
new DiffusionReaction::Jacobian<Scalar>(gto, matprop, mesh, geom_type));
775 this->add_vector_form(
new DiffusionReaction::Residual<Scalar>(gto, matprop, mesh, geom_type));
777 for (
unsigned int gfrom = 0; gfrom < G; gfrom++)
779 if(Ss_nnz[gto][gfrom])
781 this->add_matrix_form(
new Scattering::Jacobian<Scalar>(gto, gfrom, matprop, mesh, geom_type));
782 this->add_vector_form(
new Scattering::Residual<Scalar>(gto, gfrom, matprop, mesh, geom_type));
787 this->add_matrix_form(
new FissionYield::Jacobian<Scalar>(gto, gfrom, matprop, mesh, geom_type));
788 this->add_vector_form(
new FissionYield::Residual<Scalar>(gto, gfrom, matprop, mesh, geom_type));
794 template<
typename Scalar>
795 DefaultWeakFormFixedSource<Scalar>::DefaultWeakFormFixedSource(
const MaterialPropertyMaps& matprop, Mesh *mesh,
796 GeomType geom_type) : WeakForm<Scalar>(matprop.get_G())
798 lhs_init(matprop.get_G(), matprop, mesh, geom_type);
799 for (
unsigned int gto = 0; gto < matprop.get_G(); gto++)
800 this->add_vector_form(
new ExternalSources::LinearForm<Scalar>(gto, matprop, mesh, geom_type));
803 template<
typename Scalar>
804 DefaultWeakFormFixedSource<Scalar>::DefaultWeakFormFixedSource(
const MaterialPropertyMaps& matprop, Mesh *mesh,
805 Hermes2DFunction<Scalar>*f_src,
std::string src_area,
806 GeomType geom_type ) : WeakForm<Scalar>(matprop.get_G())
808 lhs_init(matprop.get_G(), matprop, mesh, geom_type);
809 for (
unsigned int gto = 0; gto < matprop.get_G(); gto++)
810 this->add_vector_form(
new WeakFormsH1::DefaultVectorFormVol<Scalar>(gto, src_area, f_src, geom_type));
813 template<
typename Scalar>
814 DefaultWeakFormFixedSource<Scalar>::DefaultWeakFormFixedSource(
const MaterialPropertyMaps& matprop, Mesh *mesh,
815 Hermes2DFunction<Scalar>*f_src,
816 Hermes::vector<std::string> src_areas,
817 GeomType geom_type ) : WeakForm<Scalar>(matprop.get_G())
819 lhs_init(matprop.get_G(), matprop, mesh, geom_type);
820 for (
unsigned int gto = 0; gto < matprop.get_G(); gto++)
821 this->add_vector_form(
new WeakFormsH1::DefaultVectorFormVol<Scalar>(gto, src_areas, f_src, geom_type));
824 template<
typename Scalar>
825 DefaultWeakFormFixedSource<Scalar>::DefaultWeakFormFixedSource(
const MaterialPropertyMaps& matprop, Mesh *mesh,
826 const Hermes::vector<Hermes2DFunction<Scalar>*>& f_src,
828 GeomType geom_type ) : WeakForm<Scalar>(matprop.get_G())
830 if(f_src.size() != matprop.get_G())
831 throw Hermes::Exceptions::Exception(E_INVALID_SIZE);
833 lhs_init(matprop.get_G(), matprop, mesh, geom_type);
834 for (
unsigned int gto = 0; gto < matprop.get_G(); gto++)
835 this->add_vector_form(
new WeakFormsH1::DefaultVectorFormVol<Scalar>(gto, src_area, f_src[gto], geom_type));
838 template<
typename Scalar>
839 DefaultWeakFormFixedSource<Scalar>::DefaultWeakFormFixedSource(
const MaterialPropertyMaps& matprop, Mesh *mesh,
840 const Hermes::vector<Hermes2DFunction<Scalar>*>& f_src,
841 Hermes::vector<std::string> src_areas,
842 GeomType geom_type ) : WeakForm<Scalar>(matprop.get_G())
844 if(f_src.size() != matprop.get_G())
845 throw Hermes::Exceptions::Exception(E_INVALID_SIZE);
847 lhs_init(matprop.get_G(), matprop, mesh, geom_type);
848 for (
unsigned int gto = 0; gto < matprop.get_G(); gto++)
849 this->add_vector_form(
new WeakFormsH1::DefaultVectorFormVol<Scalar>(gto, src_areas, f_src[gto], geom_type));
852 template<
typename Scalar>
853 DefaultWeakFormSourceIteration<Scalar>::DefaultWeakFormSourceIteration(
const MaterialPropertyMaps& matprop, Mesh *mesh,
854 Hermes::vector<MeshFunction<Scalar>*>& iterates,
855 double initial_keff_guess,
856 GeomType geom_type ) : WeakForm<Scalar>(matprop.get_G())
858 bool2 Ss_nnz = matprop.get_scattering_multigroup_structure();
860 for (
unsigned int gto = 0; gto < matprop.get_G(); gto++)
862 this->add_matrix_form(
new DiffusionReaction::Jacobian<Scalar>(gto, matprop, mesh, geom_type));
863 this->add_vector_form(
new DiffusionReaction::Residual<Scalar>(gto, matprop, mesh, geom_type));
865 for (
unsigned int gfrom = 0; gfrom < matprop.get_G(); gfrom++)
867 if(Ss_nnz[gto][gfrom])
869 this->add_matrix_form(
new Scattering::Jacobian<Scalar>(gto, gfrom, matprop, mesh, geom_type));
870 this->add_vector_form(
new Scattering::Residual<Scalar>(gto, gfrom, matprop, mesh, geom_type));
874 FissionYield::OuterIterationForm<Scalar>* keff_iteration_form =
875 new FissionYield::OuterIterationForm<Scalar>( gto, matprop, mesh, iterates, initial_keff_guess, geom_type );
876 keff_iteration_forms.push_back(keff_iteration_form);
877 this->add_vector_form(keff_iteration_form);
881 template<
typename Scalar>
882 void DefaultWeakFormSourceIteration<Scalar>::update_keff(
double new_keff)
889 for(
int i = 0; i < keff_iteration_forms.size(); i++)
890 keff_iteration_forms[i]->update_keff(new_keff);
895 namespace SupportClasses
897 void SourceFilter::filter_fn(
int n, Hermes::vector<double*> values,
double* result)
899 for (
int i = 0; i < n; i++)
902 for (
unsigned int j = 0; j < values.size(); j++)
903 result[i] += nu[j] * Sigma_f[j] * values.at(j)[i];
908 namespace Monoenergetic
912 template class HERMES_API DefaultWeakFormFixedSource<double>;
913 template class HERMES_API DefaultWeakFormFixedSource<std::complex<double> >;
918 namespace ElementaryForms
922 template class HERMES_API FissionYield::Jacobian<double>;
923 template class HERMES_API FissionYield::Jacobian<std::complex<double> >;
924 template class HERMES_API FissionYield::OuterIterationForm<double>;
925 template class HERMES_API FissionYield::OuterIterationForm<std::complex<double> >;
926 template class HERMES_API FissionYield::Residual<double>;
927 template class HERMES_API FissionYield::Residual<std::complex<double> >;
929 template class HERMES_API VacuumBoundaryCondition::Jacobian<double>;
930 template class HERMES_API VacuumBoundaryCondition::Jacobian<std::complex<double> >;
931 template class HERMES_API VacuumBoundaryCondition::Residual<double>;
932 template class HERMES_API VacuumBoundaryCondition::Residual<std::complex<double> >;
934 template double VacuumBoundaryCondition::Jacobian<double>::matrix_form<double,
double>(
int n,
double *wt, Func<double> *u_ext[], Func<double> *u,
935 Func<double> *v, Geom<double> *e, Func<double> **ext)
const;
937 template double VacuumBoundaryCondition::Residual<double>::vector_form<double,
double>(
int n,
double *wt, Func<double> *u_ext[],
938 Func<double> *v, Geom<double> *e, Func<double> **ext)
const;
940 template double DiffusionReaction::Jacobian<double>::matrix_form<double,
double>(
int n,
double *wt, Func<double> *u_ext[], Func<double> *u,
941 Func<double> *v, Geom<double> *e, Func<double> **ext)
const;
943 template double DiffusionReaction::Residual<double>::vector_form<double,
double>(
int n,
double *wt, Func<double> *u_ext[],
944 Func<double> *v, Geom<double> *e, Func<double> **ext)
const;
946 template double FissionYield::Jacobian<double>::matrix_form<double,
double>(
int n,
double *wt, Func<double> *u_ext[], Func<double> *u,
947 Func<double> *v, Geom<double> *e, Func<double> **ext)
const;
949 template double FissionYield::OuterIterationForm<double>::vector_form<double,
double>(
int n,
double *wt, Func<double> *u_ext[],
950 Func<double> *v, Geom<double> *e, Func<double> **ext )
const;
952 template double FissionYield::Residual<double>::vector_form<double,
double>(
int n,
double *wt, Func<double> *u_ext[],
953 Func<double> *v, Geom<double> *e, Func<double> **ext )
const;
955 template double Scattering::Jacobian<double>::matrix_form<double,
double>(
int n,
double *wt, Func<double> *u_ext[], Func<double> *u,
956 Func<double> *v, Geom<double> *e, Func<double> **ext )
const;
958 template double Scattering::Residual<double>::vector_form<double,
double>(
int n,
double *wt, Func<double> *u_ext[],
959 Func<double> *v, Geom<double> *e, Func<double> **ext )
const;
961 template double ExternalSources::LinearForm<double>::vector_form<double,
double>(
int n,
double *wt, Func<double> *u_ext[],
962 Func<double> *v, Geom<double> *e, Func<double> **ext)
const;
966 namespace CompleteWeakForms
970 template class HERMES_API DefaultWeakFormFixedSource<double>;
971 template class HERMES_API DefaultWeakFormFixedSource<std::complex<double> >;
973 template class HERMES_API DefaultWeakFormSourceIteration<double>;
974 template class HERMES_API DefaultWeakFormSourceIteration<std::complex<double> >;