Hermes2D  3.0
weakforms_maxwell.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 "weakform_library/weakforms_maxwell.h"
17 #include "weakform_library/integrals_h1.h"
18 
19 namespace Hermes
20 {
21  namespace Hermes2D
22  {
23  namespace WeakFormsMaxwell
24  {
25  template<typename Scalar>
26  DefaultJacobianMagnetostatics<Scalar>::DefaultJacobianMagnetostatics(int i, int j, std::string area, Scalar const_coeff,
27  CubicSpline* c_spline,
28  SymFlag sym,
29  GeomType gt,
30  int order_increase)
31  : MatrixFormVol<Scalar>(i, j), idx_j(j), const_coeff(const_coeff),
32  spline_coeff(c_spline), gt(gt),
33  order_increase(order_increase)
34  {
35  this->set_area(area);
36  this->setSymFlag(sym);
37 
38  // If spline is nullptr, initialize it to be constant 1.0.
39  if (c_spline == nullptr) this->spline_coeff = new CubicSpline(1.0);
40  }
41  template<typename Scalar>
42  DefaultJacobianMagnetostatics<Scalar>::DefaultJacobianMagnetostatics(int i, int j, std::vector<std::string> areas,
43  Scalar const_coeff,
44  CubicSpline* c_spline,
45  SymFlag sym,
46  GeomType gt,
47  int order_increase)
48  : MatrixFormVol<Scalar>(i, j), idx_j(j), const_coeff(const_coeff),
49  spline_coeff(c_spline), gt(gt),
50  order_increase(order_increase)
51  {
52  this->set_areas(areas);
53  this->setSymFlag(sym);
54 
55  // If spline is nullptr, initialize it to be constant 1.0.
56  if (c_spline == nullptr) this->spline_coeff = new CubicSpline(1.0);
57  }
58 
59  template<typename Scalar>
60  Scalar DefaultJacobianMagnetostatics<Scalar>::value(int n, double *wt, Func<Scalar> *u_ext[], Func<double> *u,
61  Func<double> *v, GeomVol<double> *e, Func<Scalar> **ext) const
62  {
63  Scalar planar_part = 0;
64  Scalar axisym_part = 0;
65  for (int i = 0; i < n; i++)
66  {
67  Scalar B_i = sqrt(sqr(u_ext[idx_j]->dx[i]) + sqr(u_ext[idx_j]->dy[i]));
68  if (std::abs(B_i) > Hermes::HermesSqrtEpsilon)
69  {
70  planar_part += wt[i] * const_coeff*spline_coeff->derivative(B_i) / B_i
71  * (u_ext[idx_j]->dx[i] * u->dx[i] + u_ext[idx_j]->dy[i] * u->dy[i])
72  * (u_ext[idx_j]->dx[i] * v->dx[i] + u_ext[idx_j]->dy[i] * v->dy[i]);
73  if (gt == HERMES_AXISYM_X)
74  {
75  axisym_part += wt[i] * const_coeff*spline_coeff->derivative(B_i) / B_i / e->y[i]
76  * (u_ext[idx_j]->dx[i] * u->dx[i] + u_ext[idx_j]->dy[i] * u->dy[i])
77  * u_ext[idx_j]->val[i] * v->dy[i];
78  }
79  else if (gt == HERMES_AXISYM_Y)
80  {
81  axisym_part += wt[i] * const_coeff*spline_coeff->derivative(B_i) / B_i / e->x[i]
82  * (u_ext[idx_j]->dx[i] * u->dx[i] + u_ext[idx_j]->dy[i] * u->dy[i])
83  * u_ext[idx_j]->val[i] * v->dx[i];
84  }
85  }
86  planar_part += wt[i] * const_coeff*spline_coeff->value(B_i)
87  * (u->dx[i] * v->dx[i] + u->dy[i] * v->dy[i]);
88  if (gt == HERMES_AXISYM_X)
89  {
90  axisym_part += wt[i] * const_coeff*spline_coeff->value(B_i) / e->y[i]
91  * u->val[i] * v->dy[i];
92  }
93  else if (gt == HERMES_AXISYM_Y)
94  {
95  axisym_part += wt[i] * const_coeff*spline_coeff->value(B_i) / e->x[i]
96  * u->val[i] * v->dx[i];
97  }
98  }
99 
100  return planar_part + axisym_part;
101  }
102 
103  template<typename Scalar>
104  Ord DefaultJacobianMagnetostatics<Scalar>::ord(int n, double *wt, Func<Ord> *u_ext[], Func<Ord> *u, Func<Ord> *v,
105  GeomVol<Ord> *e, Func<Ord> **ext) const
106  {
107  Ord planar_part(0);
108  for (int i = 0; i < n; i++)
109  {
110  Ord B_i = sqrt(sqr(u_ext[idx_j]->dx[i]) + sqr(u_ext[idx_j]->dy[i]));
111  planar_part += wt[i] * const_coeff*spline_coeff->derivative(B_i) / B_i
112  * (u_ext[idx_j]->dx[i] * u->dx[i] + u_ext[idx_j]->dy[i] * u->dy[i])
113  * (u_ext[idx_j]->dx[i] * v->dx[i] + u_ext[idx_j]->dy[i] * v->dy[i]);
114  planar_part += wt[i] * const_coeff*spline_coeff->value(B_i)
115  * (u->dx[i] * v->dx[i] + u->dy[i] * v->dy[i]);
116  }
117 
118  // This increase is for the axisymmetric part. We are not letting the
119  // Ord class do it since it would automatically choose the highest order
120  // due to the nonpolynomial 1/r term.
121  return planar_part * Ord(order_increase);
122  }
123 
124  template<typename Scalar>
125  MatrixFormVol<Scalar>* DefaultJacobianMagnetostatics<Scalar>::clone() const
126  {
127  return new DefaultJacobianMagnetostatics(*this);
128  }
129 
130  template<typename Scalar>
131  DefaultResidualMagnetostatics<Scalar>::DefaultResidualMagnetostatics(int i, std::string area, Scalar const_coeff,
132  CubicSpline* c_spline,
133  GeomType gt,
134  int order_increase)
135  : VectorFormVol<Scalar>(i), idx_i(i), const_coeff(const_coeff), spline_coeff(c_spline),
136  gt(gt), order_increase(order_increase)
137  {
138  this->set_area(area);
139  // If spline is nullptr, initialize it to be constant 1.0.
140  if (c_spline == nullptr) this->spline_coeff = new CubicSpline(1.0);
141  }
142 
143  template<typename Scalar>
144  DefaultResidualMagnetostatics<Scalar>::DefaultResidualMagnetostatics(int i, std::vector<std::string> areas, Scalar const_coeff,
145  CubicSpline* c_spline,
146  GeomType gt, int order_increase)
147  : VectorFormVol<Scalar>(i), idx_i(i), const_coeff(const_coeff), spline_coeff(c_spline), gt(gt),
148  order_increase(order_increase)
149  {
150  this->set_areas(areas);
151 
152  // If spline is nullptr, initialize it to be constant 1.0.
153  if (c_spline == nullptr) this->spline_coeff = new CubicSpline(1.0);
154  }
155 
156  template<typename Scalar>
157  Scalar DefaultResidualMagnetostatics<Scalar>::value(int n, double *wt, Func<Scalar> *u_ext[], Func<double> *v,
158  GeomVol<double> *e, Func<Scalar> **ext) const
159  {
160  Scalar planar_part = 0;
161  Scalar axisym_part = 0;
162  for (int i = 0; i < n; i++)
163  {
164  Scalar B_i = sqrt(sqr(u_ext[idx_i]->dx[i]) + sqr(u_ext[idx_i]->dy[i]));
165  planar_part += wt[i] * const_coeff*spline_coeff->value(B_i) *
166  (u_ext[idx_i]->dx[i] * v->dx[i] + u_ext[idx_i]->dy[i] * v->dy[i]);
167  if (gt == HERMES_AXISYM_X) axisym_part += wt[i] * const_coeff*spline_coeff->value(B_i) / e->y[i]
168  * u_ext[idx_i]->val[i] * v->dy[i];
169  else if (gt == HERMES_AXISYM_Y) axisym_part += wt[i] * const_coeff*spline_coeff->value(B_i) / e->x[i]
170  * u_ext[idx_i]->val[i] * v->dx[i];
171  }
172  return planar_part + axisym_part;
173  }
174 
175  template<typename Scalar>
176  Ord DefaultResidualMagnetostatics<Scalar>::ord(int n, double *wt, Func<Ord> *u_ext[], Func<Ord> *v,
177  GeomVol<Ord> *e, Func<Ord> **ext) const
178  {
179  Ord planar_part(0);
180  for (int i = 0; i < n; i++)
181  {
182  Ord B_i = sqrt(sqr(u_ext[idx_i]->dx[i]) + sqr(u_ext[idx_i]->dy[i]));
183  planar_part += wt[i] * const_coeff*spline_coeff->value(B_i) *
184  (u_ext[idx_i]->dx[i] * v->dx[i] + u_ext[idx_i]->dy[i] * v->dy[i]);
185  }
186  return planar_part * Ord(order_increase);
187  }
188 
189  // This is to make the form usable in rk_time_step_newton().
190  template<typename Scalar>
191  VectorFormVol<Scalar>* DefaultResidualMagnetostatics<Scalar>::clone() const
192  {
193  return new DefaultResidualMagnetostatics(*this);
194  }
195 
196  template class HERMES_API DefaultJacobianMagnetostatics < double > ;
197  template class HERMES_API DefaultResidualMagnetostatics < double > ;
198  }
199  }
200 }
Definition: adapt.h:24
void set_area(std::string area)
Definition: weakform.cpp:326
SymFlag
Bilinear form symmetry flag, see WeakForm::add_matrix_form.
Definition: global.h:156
::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.
Definition: global.h:148