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