18 #ifndef __H2D_OGPROJECTION_NOX_H
19 #define __H2D_OGPROJECTION_NOX_H
21 #include "../function/solution.h"
23 #include "../weakform/weakform.h"
26 #if(defined HAVE_NOX && defined HAVE_EPETRA && defined HAVE_TEUCHOS)
28 #ifdef _POSIX_C_SOURCE
29 # undef _POSIX_C_SOURCE // pyconfig.h included by NOX_Epetra defines it
32 # undef _XOPEN_SOURCE // pyconfig.h included by NOX_Epetra defines it
34 #include <NOX_Epetra.H>
41 template<
typename Scalar>
42 class HERMES_API OGProjectionNOX :
public Hermes::Mixins::Loggable
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);
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);
78 void project_global(
const Space<Scalar>* space,
79 Solution<Scalar>* source_sln, Solution<Scalar>* target_sln,
81 double newton_tol = 1e-6,
int newton_max_iter = 10);
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);
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);
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);
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);
108 class ProjectionMatrixFormVol :
public MatrixFormVol<Scalar>
111 ProjectionMatrixFormVol(
int i,
int j,
ProjNormType projNormType) : MatrixFormVol<Scalar>(i, j)
113 this->projNormType = projNormType;
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
119 switch (projNormType)
122 return l2_projection_biform<double, Scalar>(n, wt, u_ext, u, v, e, ext);
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);
132 throw Hermes::Exceptions::Exception(
"Unknown projection type");
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
140 switch (projNormType)
143 return l2_projection_biform<Hermes::Ord, Hermes::Ord>(n, wt, u_ext, u, v, e, ext);
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);
153 throw Hermes::Exceptions::Exception(
"Unknown projection type");
154 return Hermes::Ord();
158 MatrixFormVol<Scalar>* clone()
160 return new ProjectionMatrixFormVol(*
this);
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)
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]);
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)
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]);
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)
190 SolFunctionDomain result = SolFunctionDomain(0);
191 for (
int i = 0; i < n; i++)
192 result += wt[i] * (u->val[i] * v->val[i]);
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)
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]));
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)
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]));
222 class ProjectionVectorFormVol :
public VectorFormVol<Scalar>
225 ProjectionVectorFormVol(
int i,
ProjNormType projNormType) : VectorFormVol<Scalar>(i)
227 this->projNormType = projNormType;
230 Scalar value(
int n,
double *wt, Func<Scalar> *u_ext[], Func<double> *v,
231 Geom<double> *e, Func<Scalar> **ext)
const
233 switch (projNormType)
236 return l2_projection_residual<double, Scalar>(n, wt, u_ext, v, e, ext);
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);
246 throw Hermes::Exceptions::Exception(
"Unknown projection type");
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
254 switch (projNormType)
257 return l2_projection_residual<Hermes::Ord, Hermes::Ord>(n, wt, u_ext, v, e, ext);
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);
267 throw Hermes::Exceptions::Exception(
"Unknown projection type");
268 return Hermes::Ord();
272 VectorFormVol<Scalar>* clone()
274 return new ProjectionVectorFormVol(*
this);
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
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]);
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
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]);
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
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];
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
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]));
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
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]));