Hermes2D  2.0
ogprojection_nox.h
Go to the documentation of this file.
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/>.
18 #ifndef __H2D_OGPROJECTION_NOX_H
19 #define __H2D_OGPROJECTION_NOX_H
20 
21 #include "../function/solution.h"
22 #include "../forms.h"
23 #include "../weakform/weakform.h"
24 
25 //#include "epetra.h"
26 #if(defined HAVE_NOX && defined HAVE_EPETRA && defined HAVE_TEUCHOS)
27 #include <NOX.H>
28 #ifdef _POSIX_C_SOURCE
29 # undef _POSIX_C_SOURCE // pyconfig.h included by NOX_Epetra defines it
30 #endif
31 #ifdef _XOPEN_SOURCE
32 # undef _XOPEN_SOURCE // pyconfig.h included by NOX_Epetra defines it
33 #endif
34 #include <NOX_Epetra.H>
35 
36 namespace Hermes
37 {
38  namespace Hermes2D
39  {
41  template<typename Scalar>
42  class HERMES_API OGProjectionNOX : public Hermes::Mixins::Loggable
43  {
44  public:
45  OGProjectionNOX();
46 
49  void project_global(const Space<Scalar>* space,
50  MatrixFormVol<Scalar>* custom_projection_jacobian,
51  VectorFormVol<Scalar>* custom_projection_residual,
52  Scalar* target_vec, double newton_tol = 1e-6, int newton_max_iter = 10);
53 
73  void project_global(const Space<Scalar>* space, MeshFunction<Scalar>* source_meshfn,
74  Scalar* target_vec, ProjNormType proj_norm = HERMES_UNSET_NORM,
75  double newton_tol = 1e-6, int newton_max_iter = 10);
76 
78  void project_global(const Space<Scalar>* space,
79  Solution<Scalar>* source_sln, Solution<Scalar>* target_sln,
80  ProjNormType proj_norm = HERMES_UNSET_NORM,
81  double newton_tol = 1e-6, int newton_max_iter = 10);
82 
84  void project_global(Hermes::vector<const Space<Scalar>*> spaces, Hermes::vector<MeshFunction<Scalar>*> source_meshfns,
85  Scalar* target_vec, Hermes::vector<ProjNormType> proj_norms = Hermes::vector<ProjNormType>(),
86  double newton_tol = 1e-6, int newton_max_iter = 10);
87 
89  void project_global(Hermes::vector<const Space<Scalar>*> spaces, Hermes::vector<Solution<Scalar>*> source_slns,
90  Scalar* target_vec, Hermes::vector<ProjNormType> proj_norms = Hermes::vector<ProjNormType>(),
91  double newton_tol = 1e-6, int newton_max_iter = 10);
92 
93  void project_global(Hermes::vector<const Space<Scalar>*> spaces,
94  Hermes::vector<Solution<Scalar>*> source_slns, Hermes::vector<Solution<Scalar>*> target_slns,
95  Hermes::vector<ProjNormType> proj_norms = Hermes::vector<ProjNormType>(), bool delete_old_mesh = false,
96  double newton_tol = 1e-6, int newton_max_iter = 10);
97 
98  protected:
104  void project_internal(const Space<Scalar>* space, WeakForm<Scalar>* proj_wf, Scalar* target_vec,
105  double newton_tol = 1e-6, int newton_max_iter = 10);
106 
108  class ProjectionMatrixFormVol : public MatrixFormVol<Scalar>
109  {
110  public:
111  ProjectionMatrixFormVol(int i, int j, ProjNormType projNormType) : MatrixFormVol<Scalar>(i, j)
112  {
113  this->projNormType = projNormType;
114  }
115 
116  Scalar value(int n, double *wt, Func<Scalar> *u_ext[], Func<double> *u, Func<double> *v,
117  Geom<double> *e, Func<Scalar> **ext) const
118  {
119  switch (projNormType)
120  {
121  case HERMES_L2_NORM:
122  return l2_projection_biform<double, Scalar>(n, wt, u_ext, u, v, e, ext);
123  case HERMES_H1_NORM:
124  return h1_projection_biform<double, Scalar>(n, wt, u_ext, u, v, e, ext);
125  case HERMES_H1_SEMINORM:
126  return h1_semi_projection_biform<double, Scalar>(n, wt, u_ext, u, v, e, ext);
127  case HERMES_HCURL_NORM:
128  return hcurl_projection_biform<double, Scalar>(n, wt, u_ext, u, v, e, ext);
129  case HERMES_HDIV_NORM:
130  return hdiv_projection_biform<double, Scalar>(n, wt, u_ext, u, v, e, ext);
131  default:
132  throw Hermes::Exceptions::Exception("Unknown projection type");
133  return 0.0;
134  }
135  }
136 
137  Hermes::Ord ord(int n, double *wt, Func<Hermes::Ord> *u_ext[], Func<Hermes::Ord> *u, Func<Hermes::Ord> *v,
138  Geom<Hermes::Ord> *e, Func<Ord> **ext) const
139  {
140  switch (projNormType)
141  {
142  case HERMES_L2_NORM:
143  return l2_projection_biform<Hermes::Ord, Hermes::Ord>(n, wt, u_ext, u, v, e, ext);
144  case HERMES_H1_NORM:
145  return h1_projection_biform<Hermes::Ord, Hermes::Ord>(n, wt, u_ext, u, v, e, ext);
146  case HERMES_H1_SEMINORM:
147  return h1_semi_projection_biform<Hermes::Ord, Hermes::Ord>(n, wt, u_ext, u, v, e, ext);
148  case HERMES_HCURL_NORM:
149  return hcurl_projection_biform<Hermes::Ord, Hermes::Ord>(n, wt, u_ext, u, v, e, ext);
150  case HERMES_HDIV_NORM:
151  return hdiv_projection_biform<Hermes::Ord, Hermes::Ord>(n, wt, u_ext, u, v, e, ext);
152  default:
153  throw Hermes::Exceptions::Exception("Unknown projection type");
154  return Hermes::Ord();
155  }
156  }
157 
158  MatrixFormVol<Scalar>* clone()
159  {
160  return new ProjectionMatrixFormVol(*this);
161  }
162 
163  private:
164  ProjNormType projNormType;
165 
166  template<typename TestFunctionDomain, typename SolFunctionDomain>
167  static SolFunctionDomain h1_projection_biform(int n, double *wt, Func<SolFunctionDomain> *u_ext[], Func<TestFunctionDomain> *u,
168  Func<TestFunctionDomain> *v, Geom<TestFunctionDomain> *e, Func<SolFunctionDomain> **ext)
169  {
170  SolFunctionDomain result = SolFunctionDomain(0);
171  for (int i = 0; i < n; i++)
172  result += wt[i] * (u->val[i] * v->val[i] + u->dx[i] * v->dx[i] + u->dy[i] * v->dy[i]);
173  return result;
174  }
175 
176  template<typename TestFunctionDomain, typename SolFunctionDomain>
177  static SolFunctionDomain h1_semi_projection_biform(int n, double *wt, Func<SolFunctionDomain> *u_ext[], Func<TestFunctionDomain> *u,
178  Func<TestFunctionDomain> *v, Geom<TestFunctionDomain> *e, Func<SolFunctionDomain> **ext)
179  {
180  SolFunctionDomain result = SolFunctionDomain(0);
181  for (int i = 0; i < n; i++)
182  result += wt[i] * (u->dx[i] * v->dx[i] + u->dy[i] * v->dy[i]);
183  return result;
184  }
185 
186  template<typename TestFunctionDomain, typename SolFunctionDomain>
187  static SolFunctionDomain l2_projection_biform(int n, double *wt, Func<SolFunctionDomain> *u_ext[], Func<TestFunctionDomain> *u,
188  Func<TestFunctionDomain> *v, Geom<TestFunctionDomain> *e, Func<SolFunctionDomain> **ext)
189  {
190  SolFunctionDomain result = SolFunctionDomain(0);
191  for (int i = 0; i < n; i++)
192  result += wt[i] * (u->val[i] * v->val[i]);
193  return result;
194  }
195 
196  template<typename TestFunctionDomain, typename SolFunctionDomain>
197  static SolFunctionDomain hcurl_projection_biform(int n, double *wt, Func<SolFunctionDomain> *u_ext[], Func<TestFunctionDomain> *u,
198  Func<TestFunctionDomain> *v, Geom<TestFunctionDomain> *e, Func<SolFunctionDomain> **ext)
199  {
200  SolFunctionDomain result = SolFunctionDomain(0);
201  for (int i = 0; i < n; i++) {
202  result += wt[i] * (u->curl[i] * conj(v->curl[i]));
203  result += wt[i] * (u->val0[i] * conj(v->val0[i]) + u->val1[i] * conj(v->val1[i]));
204  }
205  return result;
206  }
207 
208  template<typename TestFunctionDomain, typename SolFunctionDomain>
209  static SolFunctionDomain hdiv_projection_biform(int n, double *wt, Func<SolFunctionDomain> *u_ext[], Func<TestFunctionDomain> *u,
210  Func<TestFunctionDomain> *v, Geom<TestFunctionDomain> *e, Func<SolFunctionDomain> **ext)
211  {
212  SolFunctionDomain result = SolFunctionDomain(0);
213  for (int i = 0; i < n; i++) {
214  result += wt[i] * (u->div[i] * conj(v->div[i]));
215  result += wt[i] * (u->val0[i] * conj(v->val0[i]) + u->val1[i] * conj(v->val1[i]));
216  }
217  return result;
218  }
219  };
220 
222  class ProjectionVectorFormVol : public VectorFormVol<Scalar>
223  {
224  public:
225  ProjectionVectorFormVol(int i, ProjNormType projNormType) : VectorFormVol<Scalar>(i)
226  {
227  this->projNormType = projNormType;
228  }
229 
230  Scalar value(int n, double *wt, Func<Scalar> *u_ext[], Func<double> *v,
231  Geom<double> *e, Func<Scalar> **ext) const
232  {
233  switch (projNormType)
234  {
235  case HERMES_L2_NORM:
236  return l2_projection_residual<double, Scalar>(n, wt, u_ext, v, e, ext);
237  case HERMES_H1_NORM:
238  return h1_projection_residual<double, Scalar>(n, wt, u_ext, v, e, ext);
239  case HERMES_H1_SEMINORM:
240  return h1_semi_projection_residual<double, Scalar>(n, wt, u_ext, v, e, ext);
241  case HERMES_HCURL_NORM:
242  return hcurl_projection_residual<double, Scalar>(n, wt, u_ext, v, e, ext);
243  case HERMES_HDIV_NORM:
244  return hdiv_projection_residual<double, Scalar>(n, wt, u_ext, v, e, ext);
245  default:
246  throw Hermes::Exceptions::Exception("Unknown projection type");
247  return 0.0;
248  }
249  }
250 
251  Hermes::Ord ord(int n, double *wt, Func<Hermes::Ord> *u_ext[], Func<Hermes::Ord> *v,
252  Geom<Hermes::Ord> *e, Func<Ord> **ext) const
253  {
254  switch (projNormType)
255  {
256  case HERMES_L2_NORM:
257  return l2_projection_residual<Hermes::Ord, Hermes::Ord>(n, wt, u_ext, v, e, ext);
258  case HERMES_H1_NORM:
259  return h1_projection_residual<Hermes::Ord, Hermes::Ord>(n, wt, u_ext, v, e, ext);
260  case HERMES_H1_SEMINORM:
261  return h1_semi_projection_residual<Hermes::Ord, Hermes::Ord>(n, wt, u_ext, v, e, ext);
262  case HERMES_HCURL_NORM:
263  return hcurl_projection_residual<Hermes::Ord, Hermes::Ord>(n, wt, u_ext, v, e, ext);
264  case HERMES_HDIV_NORM:
265  return hdiv_projection_residual<Hermes::Ord, Hermes::Ord>(n, wt, u_ext, v, e, ext);
266  default:
267  throw Hermes::Exceptions::Exception("Unknown projection type");
268  return Hermes::Ord();
269  }
270  }
271 
272  VectorFormVol<Scalar>* clone()
273  {
274  return new ProjectionVectorFormVol(*this);
275  }
276 
277  private:
278  ProjNormType projNormType;
279 
280  template<typename TestFunctionDomain, typename SolFunctionDomain>
281  SolFunctionDomain h1_projection_residual(int n, double *wt, Func<SolFunctionDomain> *u_ext[], Func<TestFunctionDomain> *v,
282  Geom<TestFunctionDomain> *e, Func<SolFunctionDomain> **ext) const
283  {
284  SolFunctionDomain result = SolFunctionDomain(0);
285  for (int i = 0; i < n; i++)
286  result += wt[i] * ((u_ext[this->i]->val[i] - ext[0]->val[i]) * v->val[i]
287  + (u_ext[this->i]->dx[i] - ext[0]->dx[i]) * v->dx[i]
288  + (u_ext[this->i]->dy[i] - ext[0]->dy[i]) * v->dy[i]);
289  return result;
290  }
291 
292  template<typename TestFunctionDomain, typename SolFunctionDomain>
293  SolFunctionDomain h1_semi_projection_residual(int n, double *wt, Func<SolFunctionDomain> *u_ext[], Func<TestFunctionDomain> *v,
294  Geom<TestFunctionDomain> *e, Func<SolFunctionDomain> **ext) const
295  {
296  SolFunctionDomain result = SolFunctionDomain(0);
297  for (int i = 0; i < n; i++)
298  result += wt[i] * ((u_ext[this->i]->dx[i] - ext[0]->dx[i]) * v->dx[i]
299  + (u_ext[this->i]->dy[i] - ext[0]->dy[i]) * v->dy[i]);
300  return result;
301  }
302 
303  template<typename TestFunctionDomain, typename SolFunctionDomain>
304  SolFunctionDomain l2_projection_residual(int n, double *wt, Func<SolFunctionDomain> *u_ext[], Func<TestFunctionDomain> *v,
305  Geom<TestFunctionDomain> *e, Func<SolFunctionDomain> **ext) const
306  {
307  SolFunctionDomain result = SolFunctionDomain(0);
308  for (int i = 0; i < n; i++)
309  result += wt[i] * (u_ext[this->i]->val[i] - ext[0]->val[i]) * v->val[i];
310  return result;
311  }
312 
313  template<typename TestFunctionDomain, typename SolFunctionDomain>
314  SolFunctionDomain hcurl_projection_residual(int n, double *wt, Func<SolFunctionDomain> *u_ext[], Func<TestFunctionDomain> *v,
315  Geom<TestFunctionDomain> *e, Func<SolFunctionDomain> **ext) const
316  {
317  SolFunctionDomain result = SolFunctionDomain(0);
318  for (int i = 0; i < n; i++) {
319  result += wt[i] * (u_ext[this->i]->curl[i] - ext[0]->curl[i]) * conj(v->curl[i]);
320  result += wt[i] * ((u_ext[this->i]->val0[i] - ext[0]->val0[i]) * conj(v->val0[i])
321  + (u_ext[this->i]->val1[i] - ext[0]->val1[i]) * conj(v->val1[i]));
322  }
323 
324  return result;
325  }
326 
327  template<typename TestFunctionDomain, typename SolFunctionDomain>
328  SolFunctionDomain hdiv_projection_residual(int n, double *wt, Func<SolFunctionDomain> *u_ext[], Func<TestFunctionDomain> *v,
329  Geom<TestFunctionDomain> *e, Func<SolFunctionDomain> **ext) const
330  {
331  SolFunctionDomain result = SolFunctionDomain(0);
332  for (int i = 0; i < n; i++) {
333  result += wt[i] * (u_ext[this->i]->div[i] - ext[0]->div[i]) * conj(v->div[i]);
334  result += wt[i] * ((u_ext[this->i]->val0[i] - ext[0]->val0[i]) * conj(v->val0[i])
335  + (u_ext[this->i]->val1[i] - ext[0]->val1[i]) * conj(v->val1[i]));
336  }
337 
338  return result;
339  }
340  };
341 
342  int ndof;
343  };
344  }
345 }
346 #endif
347 #endif