Hermes2D  2.0
forms.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 
17 
18 #include "forms.h"
19 #include "api2d.h"
20 #include <complex>
21 
22 namespace Hermes
23 {
24  namespace Hermes2D
25  {
26  template<typename T>
27  Func<T>::Func(int num_gip, int num_comps) : num_gip(num_gip), nc(num_comps)
28  {
29  val = NULL;
30  dx = NULL;
31  dy = NULL;
32  laplace = NULL;
33 
34  if(this->nc > 1)
35  {
36  val0 = val1 = NULL;
37  dx0 = dx1 = NULL;
38  dy0 = dy1 = NULL;
39  curl = NULL;
40  div = NULL;
41  }
42  }
43 
44  template<typename T>
45  void Func<T>::subtract(Func<T>* func)
46  {
47  if(num_gip != func->num_gip)
48  throw Hermes::Exceptions::Exception("Unable to subtract a function due to a different number of integration points (this: %d, other: %d)", num_gip, func->num_gip);
49  if(nc != func->nc)
50  throw Hermes::Exceptions::Exception("Unable to subtract a function due to a different number of components (this: %d, other: %d)", nc, func->nc);
51 
52  subtract(this->val, func->val);
53  subtract(this->dx, func->dx);
54  subtract(this->dy, func->dy);
55 
56 #ifdef H2D_USE_SECOND_DERIVATIVES
57  subtract(this->laplace, func->laplace);
58 #endif
59 
60  if(nc > 1)
61  {
62  subtract(this->val0, func->val0);
63  subtract(this->val1, func->val1);
64  subtract(this->dx0, func->dx0);
65  subtract(this->dx1, func->dx1);
66  subtract(this->dy0, func->dy0);
67  subtract(this->dy1, func->dy1);
68  subtract(this->curl, func->curl);
69  subtract(this->div, func->div);
70  }
71  };
72 
73  template<typename T>
74  int Func<T>::get_num_gip() const
75  {
76  return num_gip;
77  }
78 
79  template<typename T>
80  void Func<T>::subtract(T* attribute, T* other_attribute)
81  {
82  if(attribute != NULL && other_attribute != NULL)
83  {
84  for(int i = 0; i < num_gip; i++)
85  attribute[i] -= other_attribute[i];
86  }
87  }
88 
89  template<typename T>
90  void Func<T>::add(Func<T>* func)
91  {
92  if(num_gip != func->num_gip)
93  throw Hermes::Exceptions::Exception("Unable to add a function due to a different number of integration points (this: %d, other: %d)", num_gip, func->num_gip);
94  if(nc != func->nc)
95  throw Hermes::Exceptions::Exception("Unable to add a function due to a different number of components (this: %d, other: %d)", nc, func->nc);
96 
97  add(this->val, func->val);
98  add(this->dx, func->dx);
99  add(this->dy, func->dy);
100 
101 #ifdef H2D_USE_SECOND_DERIVATIVES
102  add(this->laplace, func->laplace);
103 #endif
104 
105  if(nc > 1)
106  {
107  add(this->val0, func->val0);
108  add(this->val1, func->val1);
109  add(this->dx0, func->dx0);
110  add(this->dx1, func->dx1);
111  add(this->dy0, func->dy0);
112  add(this->dy1, func->dy1);
113  add(this->curl, func->curl);
114  add(this->div, func->div);
115  }
116  };
117 
118  template<typename T>
119  void Func<T>::add(T* attribute, T* other_attribute)
120  {
121  if(attribute != NULL && other_attribute != NULL)
122  {
123  if(other_attribute == NULL)
124  throw Hermes::Exceptions::Exception("Unable to add a function expansion, the desired attribute is NULL in the other function.");
125  for(int i = 0; i < num_gip; i++)
126  attribute[i] += other_attribute[i];
127  }
128  }
129 
130  template<typename T>
131  void Func<T>::free_ord()
132  {
133  delete val;
134  delete dx;
135  val = NULL;
136  dx = NULL;
137  dy = NULL;
138 
139  if(this->nc > 1)
140  {
141  val0 = val1 = NULL;
142  dx0 = dx1 = NULL;
143  dy0 = dy1 = NULL;
144  curl = NULL;
145  div = NULL;
146  }
147 
148  laplace = NULL;
149  }
150 
151  template<typename T>
152  void Func<T>::free_fn()
153  {
154  delete [] val; val = NULL;
155  delete [] dx; dx = NULL;
156  delete [] dy; dy = NULL;
157 
158 #ifdef H2D_USE_SECOND_DERIVATIVES
159  delete [] laplace;
160  laplace = NULL;
161 #endif
162 
163  if(this->nc > 1)
164  {
165  delete [] val0; delete [] val1; val0 = val1 = NULL;
166  delete [] dx0; delete [] dx1; dx0 = dx1 = NULL;
167  delete [] dy0; delete [] dy1; dy0 = dy1 = NULL;
168  delete [] curl; curl = NULL;
169  delete [] div; div = NULL;
170  }
171  }
172 
173  template<typename T>
174  DiscontinuousFunc<T>::DiscontinuousFunc(Func<T>* fn, bool support_on_neighbor, bool reverse) :
175  Func<T>(fn->num_gip, fn->nc),
176  fn_central(NULL),
177  fn_neighbor(NULL),
178  reverse_neighbor_side(reverse)
179  {
180  if(fn == NULL)
181  throw Hermes::Exceptions::Exception("Invalid arguments to DiscontinuousFunc constructor.");
182  if(support_on_neighbor)
183  fn_neighbor = fn;
184  else
185  fn_central = fn;
186  }
187 
188  template<typename T>
189  DiscontinuousFunc<T>::DiscontinuousFunc(Func<T>* fn_c, Func<T>* fn_n, bool reverse) :
190  Func<T>(fn_c->num_gip, fn_c->nc),
191  fn_central(fn_c),
192  fn_neighbor(fn_n),
193  reverse_neighbor_side(reverse)
194  {
195  if(fn_c == NULL)
196  throw Hermes::Exceptions::NullException(0);
197  if(fn_n == NULL)
198  throw Hermes::Exceptions::NullException(1);
199  if(fn_c->num_gip != fn_n->num_gip || fn_c->nc != fn_n->nc)
200  throw Hermes::Exceptions::Exception("DiscontinuousFunc must be formed by two Func's with same number of integration points and components.");
201  }
202 
203  // Explicit template specializations are needed here, general template<T> T DiscontinuousFunc<T>::zero = T(0) doesn't work.
204  template<> Hermes::Ord DiscontinuousFunc<Hermes::Ord>::zero = Hermes::Ord(0);
205  template<> double DiscontinuousFunc<double>::zero = 0.0;
206  template<> std::complex<double> DiscontinuousFunc<std::complex<double> >::zero = std::complex<double>(0);
207 
208  template<typename T>
209  T& DiscontinuousFunc<T>::get_val_central(int k) const { return (fn_central != NULL) ? fn_central->val[k] : zero; }
210 
211  template<typename T>
212  T& DiscontinuousFunc<T>::get_val_neighbor(int k) const { return (fn_neighbor != NULL) ? fn_neighbor->val[ reverse_neighbor_side ? fn_neighbor->num_gip-k-1 : k ] : zero; }
213 
214  template<typename T>
215  T& DiscontinuousFunc<T>::get_dx_central(int k) const { return (fn_central != NULL) ? fn_central->dx[k] : zero; }
216 
217  template<typename T>
218  T& DiscontinuousFunc<T>::get_dx_neighbor(int k) const { return (fn_neighbor != NULL) ? fn_neighbor->dx[ reverse_neighbor_side ? fn_neighbor->num_gip-k-1 : k ] : zero; }
219 
220  template<typename T>
221  T& DiscontinuousFunc<T>::get_dy_central(int k) const { return (fn_central != NULL) ? fn_central->dy[k] : zero; }
222 
223  template<typename T>
224  T& DiscontinuousFunc<T>::get_dy_neighbor(int k) const { return (fn_neighbor != NULL) ? fn_neighbor->dy[ reverse_neighbor_side ? fn_neighbor->num_gip-k-1 : k ] : zero; }
225 
226  template<typename T>
227  T& DiscontinuousFunc<T>::get_laplace_central(int k) { return (fn_central != NULL) ? fn_central->laplace[k] : zero; }
228 
229  template<typename T>
230  T& DiscontinuousFunc<T>::get_laplace_neighbor(int k) { return (fn_neighbor != NULL) ? fn_neighbor->laplace[ reverse_neighbor_side ? fn_neighbor->num_gip-k-1 : k ] : zero; }
231 
232  template<typename T>
233  void DiscontinuousFunc<T>::subtract(const DiscontinuousFunc<T>& func)
234  {
235  if(fn_central != NULL && func.fn_central != NULL)
236  fn_central->subtract(func.fn_central);
237  if(fn_neighbor != NULL && func.fn_neighbor != NULL)
238  fn_neighbor->subtract(func.fn_neighbor);
239  }
240 
241  template<typename T>
242  void DiscontinuousFunc<T>::free_fn()
243  {
244  if(fn_central != NULL)
245  {
246  fn_central->free_fn();
247  delete fn_central;
248  fn_central = NULL;
249  }
250  if(fn_neighbor != NULL)
251  {
252  fn_neighbor->free_fn();
253  delete fn_neighbor;
254  fn_neighbor = NULL;
255  }
256  }
257 
258  template<typename T>
259  void DiscontinuousFunc<T>::free_ord()
260  {
261  if(fn_central != NULL)
262  {
263  fn_central->free_ord();
264  delete fn_central;
265  fn_central = NULL;
266  }
267  if(fn_neighbor != NULL)
268  {
269  fn_neighbor->free_ord();
270  delete fn_neighbor;
271  fn_neighbor = NULL;
272  }
273  }
274 
275  template<typename T>
276  Geom<T>::Geom()
277  {
278  elem_marker = -1;
279  edge_marker = -1;
280  id = 0;
281  isurf = 4;
282  x = y = NULL;
283  nx = ny = NULL;
284  tx = ty = NULL;
285  }
286 
287  template<typename T>
288  void Geom<T>::free()
289  {
290  delete [] x; delete [] y;
291  delete [] tx; delete [] ty;
292  delete [] nx; delete [] ny;
293  }
294 
295  template<typename T>
296  InterfaceGeom<T>::InterfaceGeom(Geom<T>* geom, int n_marker, int n_id, T n_diam) :
297  Geom<T>(), neighb_marker(n_marker), neighb_id(n_id), neighb_diam(n_diam)
298  {
299  // Let this class expose the standard Geom interface.
300  this->edge_marker = geom->edge_marker;
301  this->elem_marker = geom->elem_marker;
302  this->id = geom->id;
303  this->isurf = geom->isurf;
304  this->diam = geom->diam;
305  this->area = geom->area;
306  this->x = geom->x;
307  this->y = geom->y;
308  this->tx = geom->tx;
309  this->ty = geom->ty;
310  this->nx = geom->nx;
311  this->ny = geom->ny;
312  this->orientation = geom->orientation;
313  this->wrapped_geom = geom;
314  }
315 
316  template<typename T>
317  void InterfaceGeom<T>::free()
318  {
319  wrapped_geom->free();
320  delete wrapped_geom;
321  }
322 
323  template<typename T>
324  void InterfaceGeom<T>::free_ord()
325  {
326  delete wrapped_geom;
327  }
328 
329  template<typename T>
330  int InterfaceGeom<T>::get_neighbor_marker() const
331  {
332  return neighb_marker;
333  }
334 
335  template<typename T>
336  int InterfaceGeom<T>::get_neighbor_id() const
337  {
338  return neighb_id;
339  }
340 
341  template<typename T>
342  T InterfaceGeom<T>::get_neighbor_diam() const
343  {
344  return neighb_diam;
345  }
346 
347  Geom<Hermes::Ord>* init_geom_ord()
348  {
349  Geom<Hermes::Ord>* e = new Geom<Hermes::Ord>;
350  static Hermes::Ord x[] = { Hermes::Ord(1) };
351  static Hermes::Ord y[] = { Hermes::Ord(1) };
352 
353  static Hermes::Ord nx[] = { Hermes::Ord(1) };
354  static Hermes::Ord ny[] = { Hermes::Ord(1) };
355 
356  static Hermes::Ord tx[] = { Hermes::Ord(1) };
357  static Hermes::Ord ty[] = { Hermes::Ord(1) };
358 
359  static Hermes::Ord diam = Hermes::Ord(1);
360 
361  e->x = x; e->y = y;
362  e->nx = nx; e->ny = ny;
363  e->tx = tx; e->ty = ty;
364  e->diam = diam;
365 
366  return e;
367  }
368 
369  Geom<double>* init_geom_vol(RefMap *rm, const int order)
370  {
371  Geom<double>* e = new Geom<double>;
372  e->diam = rm->get_active_element()->get_diameter();
373  e->area = rm->get_active_element()->get_area();
374  e->id = rm->get_active_element()->id;
375  e->elem_marker = rm->get_active_element()->marker;
376  Quad2D* quad = rm->get_quad_2d();
377  int np = quad->get_num_points(order, rm->get_active_element()->get_mode());
378  e->x = new double[np];
379  e->y = new double[np];
380  double* x = rm->get_phys_x(order);
381  double* y = rm->get_phys_y(order);
382  for (int i = 0; i < np; i++)
383  {
384  e->x[i] = x[i];
385  e->y[i] = y[i];
386  }
387  return e;
388  }
389 
390  Geom<double>* init_geom_surf(RefMap *rm, int isurf, int marker, const int order, double3*& tan)
391  {
392  Geom<double>* e = new Geom<double>;
393  e->edge_marker = marker;
394  e->elem_marker = rm->get_active_element()->marker;
395  e->diam = rm->get_active_element()->get_diameter();
396  e->area = rm->get_active_element()->get_area();
397  e->id = rm->get_active_element()->id;
398  e->isurf = isurf;
399 
400  tan = rm->get_tangent(isurf, order);
401  double* x = rm->get_phys_x(order);
402  double* y = rm->get_phys_y(order);
403 
404  Quad2D* quad = rm->get_quad_2d();
405  int np = quad->get_num_points(order, rm->get_active_element()->get_mode());
406  e->x = new double[np];
407  e->y = new double[np];
408  e->tx = new double[np];
409  e->ty = new double[np];
410  e->nx = new double[np];
411  e->ny = new double[np];
412  for (int i = 0; i < np; i++)
413  {
414  e->x[i] = x[i];
415  e->y[i] = y[i];
416  e->tx[i] = tan[i][0]; e->ty[i] = tan[i][1];
417  e->nx[i] = tan[i][1]; e->ny[i] = - tan[i][0];
418  }
419  e->orientation = rm->get_active_element()->get_edge_orientation(isurf);
420  return e;
421  }
422 
423  Func<Hermes::Ord>* init_fn_ord(const int order)
424  {
425  Hermes::Ord *d = new Hermes::Ord(order);
426  Hermes::Ord *d1 = new Hermes::Ord(order > 1 ? order - 1 : order);
427 
428  Func<Hermes::Ord>* f = new Func<Hermes::Ord>(1, 2);
429  f->val = d;
430  f->dx = d1;
431  f->dy = d1;
432 #ifdef H2D_USE_SECOND_DERIVATIVES
433  f->laplace = d;
434 #endif
435  f->val0 = f->val1 = d;
436  f->dx0 = f->dx1 = d1;
437  f->dy0 = f->dy1 = d1;
438  f->curl = d1;
439  f->div = d1;
440  return f;
441  }
442 
443  Func<double>* init_fn(PrecalcShapeset *fu, RefMap *rm, const int order)
444  {
445  int nc = fu->get_num_components();
446  SpaceType space_type = fu->get_space_type();
447  Quad2D* quad = fu->get_quad_2d();
448 
449 #ifdef H2D_USE_SECOND_DERIVATIVES
450  if(space_type == HERMES_H1_SPACE || space_type == HERMES_L2_SPACE)
451  fu->set_quad_order(order, H2D_FN_ALL);
452  else
453  fu->set_quad_order(order);
454 #else
455  fu->set_quad_order(order);
456 #endif
457 
458  double3* pt = quad->get_points(order, rm->get_active_element()->get_mode());
459  int np = quad->get_num_points(order, rm->get_active_element()->get_mode());
460  Func<double>* u = new Func<double>(np, nc);
461 
462  // H1 & L2 space.
463  if(space_type == HERMES_H1_SPACE || space_type == HERMES_L2_SPACE)
464  {
465  u->val = new double[np];
466  u->dx = new double[np];
467  u->dy = new double[np];
468 
469 #ifdef H2D_USE_SECOND_DERIVATIVES
470  u->laplace = new double[np];
471 #endif
472 
473  double *fn = fu->get_fn_values();
474  double *dx = fu->get_dx_values();
475  double *dy = fu->get_dy_values();
476 
477  double *dxx;
478  double *dxy;
479  double *dyy;
480 
481 #ifdef H2D_USE_SECOND_DERIVATIVES
482  dxx = fu->get_dxx_values();
483  dxy = fu->get_dxy_values();
484  dyy = fu->get_dyy_values();
485 #endif
486 
487  double2x2 *m;
488  if(rm->is_jacobian_const())
489  {
490  m = new double2x2[np];
491  double2x2 const_inv_ref_map;
492 
493  const_inv_ref_map[0][0] = rm->get_const_inv_ref_map()[0][0][0];
494  const_inv_ref_map[0][1] = rm->get_const_inv_ref_map()[0][0][1];
495  const_inv_ref_map[1][0] = rm->get_const_inv_ref_map()[0][1][0];
496  const_inv_ref_map[1][1] = rm->get_const_inv_ref_map()[0][1][1];
497 
498  for(int i = 0; i < np; i++)
499  {
500  m[i][0][0] = const_inv_ref_map[0][0];
501  m[i][0][1] = const_inv_ref_map[0][1];
502  m[i][1][0] = const_inv_ref_map[1][0];
503  m[i][1][1] = const_inv_ref_map[1][1];
504  }
505  }
506  else
507  m = rm->get_inv_ref_map(order);
508 
509 #ifdef H2D_USE_SECOND_DERIVATIVES
510  double3x2 *mm;
511  mm = rm->get_second_ref_map(order);
512  for (int i = 0; i < np; i++, m++, mm++)
513  {
514  u->val[i] = fn[i];
515  u->dx[i] = (dx[i] * (*m)[0][0] + dy[i] * (*m)[0][1]);
516  u->dy[i] = (dx[i] * (*m)[1][0] + dy[i] * (*m)[1][1]);
517 
518  double axx = (Hermes::sqr((*m)[0][0]) + Hermes::sqr((*m)[1][0]));
519  double ayy = (Hermes::sqr((*m)[0][1]) + Hermes::sqr((*m)[1][1]));
520  double axy = 2.0 * ((*m)[0][0]*(*m)[0][1] + (*m)[1][0]*(*m)[1][1]);
521  double ax = (*mm)[0][0] + (*mm)[2][0];
522  double ay = (*mm)[0][1] + (*mm)[2][1];
523  u->laplace[i] = ( dx[i] * ax + dy[i] * ay + dxx[i] * axx + dxy[i] * axy + dyy[i] * ayy );
524  }
525 #else
526  for (int i = 0; i < np; i++, m++)
527  {
528  u->val[i] = fn[i];
529  u->dx[i] = (dx[i] * (*m)[0][0] + dy[i] * (*m)[0][1]);
530  u->dy[i] = (dx[i] * (*m)[1][0] + dy[i] * (*m)[1][1]);
531  }
532 #endif
533 
534  m -= np;
535  if(rm->is_jacobian_const())
536  delete [] m;
537  }
538  // Hcurl space.
539  else if(space_type == HERMES_HCURL_SPACE)
540  {
541  u->val0 = new double[np];
542  u->val1 = new double[np];
543  u->curl = new double[np];
544 
545  double *fn0 = fu->get_fn_values(0);
546  double *fn1 = fu->get_fn_values(1);
547  double *dx1 = fu->get_dx_values(1);
548  double *dy0 = fu->get_dy_values(0);
549  double2x2 *m;
550  if(rm->is_jacobian_const())
551  {
552  m = new double2x2[np];
553  double2x2 const_inv_ref_map;
554 
555  const_inv_ref_map[0][0] = rm->get_const_inv_ref_map()[0][0][0];
556  const_inv_ref_map[0][1] = rm->get_const_inv_ref_map()[0][0][1];
557  const_inv_ref_map[1][0] = rm->get_const_inv_ref_map()[0][1][0];
558  const_inv_ref_map[1][1] = rm->get_const_inv_ref_map()[0][1][1];
559 
560  for(int i = 0; i < np; i++)
561  {
562  m[i][0][0] = const_inv_ref_map[0][0];
563  m[i][0][1] = const_inv_ref_map[0][1];
564  m[i][1][0] = const_inv_ref_map[1][0];
565  m[i][1][1] = const_inv_ref_map[1][1];
566  }
567  }
568  else
569  m = rm->get_inv_ref_map(order);
570  for (int i = 0; i < np; i++, m++)
571  {
572  u->val0[i] = (fn0[i] * (*m)[0][0] + fn1[i] * (*m)[0][1]);
573  u->val1[i] = (fn0[i] * (*m)[1][0] + fn1[i] * (*m)[1][1]);
574  u->curl[i] = ((*m)[0][0] * (*m)[1][1] - (*m)[1][0] * (*m)[0][1]) * (dx1[i] - dy0[i]);
575  }
576 
577  m -= np;
578  if(rm->is_jacobian_const())
579  delete [] m;
580  }
581  // Hdiv space.
582  else if(space_type == HERMES_HDIV_SPACE)
583  {
584  u->val0 = new double[np];
585  u->val1 = new double[np];
586 
587  double *fn0 = fu->get_fn_values(0);
588  double *fn1 = fu->get_fn_values(1);
589  double *dx0 = fu->get_dx_values(0);
590  double *dy1 = fu->get_dy_values(1);
591  double2x2 *m;
592  if(rm->is_jacobian_const())
593  {
594  m = new double2x2[np];
595  double2x2 const_inv_ref_map;
596 
597  const_inv_ref_map[0][0] = rm->get_const_inv_ref_map()[0][0][0];
598  const_inv_ref_map[0][1] = rm->get_const_inv_ref_map()[0][0][1];
599  const_inv_ref_map[1][0] = rm->get_const_inv_ref_map()[0][1][0];
600  const_inv_ref_map[1][1] = rm->get_const_inv_ref_map()[0][1][1];
601 
602  for(int i = 0; i < np; i++)
603  {
604  m[i][0][0] = const_inv_ref_map[0][0];
605  m[i][0][1] = const_inv_ref_map[0][1];
606  m[i][1][0] = const_inv_ref_map[1][0];
607  m[i][1][1] = const_inv_ref_map[1][1];
608  }
609  }
610  else
611  m = rm->get_inv_ref_map(order);
612  for (int i = 0; i < np; i++, m++)
613  {
614  u->val0[i] = ( fn0[i] * (*m)[1][1] - fn1[i] * (*m)[1][0]);
615  u->val1[i] = (- fn0[i] * (*m)[0][1] + fn1[i] * (*m)[0][0]);
616  u->div[i] = ((*m)[0][0] * (*m)[1][1] - (*m)[1][0] * (*m)[0][1]) * (dx0[i] + dy1[i]);
617  }
618  m -= np;
619  if(rm->is_jacobian_const())
620  delete [] m;
621  }
622  else
623  throw Hermes::Exceptions::Exception("Wrong space type - space has to be either H1, Hcurl, Hdiv or L2");
624 
625  return u;
626  }
627 
628  template<typename Scalar>
629  Func<Scalar>* init_fn(MeshFunction<Scalar>*fu, const int order)
630  {
631  // Sanity checks.
632  if(fu == NULL) throw Hermes::Exceptions::Exception("NULL MeshFunction in Func<Scalar>*::init_fn().");
633  if(fu->get_mesh() == NULL) throw Hermes::Exceptions::Exception("Uninitialized MeshFunction used.");
634 
635  int nc = fu->get_num_components();
636  Quad2D* quad = fu->get_quad_2d();
637  fu->set_quad_order(order);
638  double3* pt = quad->get_points(order, fu->get_active_element()->get_mode());
639  int np = quad->get_num_points(order, fu->get_active_element()->get_mode());
640  Func<Scalar>* u = new Func<Scalar>(np, nc);
641 
642  if(u->nc == 1)
643  {
644  u->val = new Scalar[np];
645  u->dx = new Scalar[np];
646  u->dy = new Scalar[np];
647  memcpy(u->val, fu->get_fn_values(), np * sizeof(Scalar));
648  memcpy(u->dx, fu->get_dx_values(), np * sizeof(Scalar));
649  memcpy(u->dy, fu->get_dy_values(), np * sizeof(Scalar));
650  }
651  else if(u->nc == 2)
652  {
653  u->val0 = new Scalar[np];
654  u->val1 = new Scalar[np];
655  u->curl = new Scalar[np];
656  u->div = new Scalar[np];
657 
658  memcpy(u->val0, fu->get_fn_values(0), np * sizeof(Scalar));
659  memcpy(u->val1, fu->get_fn_values(1), np * sizeof(Scalar));
660 
661  Scalar *dx1 = fu->get_dx_values(1);
662  Scalar *dy0 = fu->get_dy_values(0);
663  for (int i = 0; i < np; i++) u->curl[i] = dx1[i] - dy0[i];
664 
665  Scalar *dx0 = fu->get_dx_values(0);
666  Scalar *dy1 = fu->get_dy_values(1);
667  for (int i = 0; i < np; i++) u->div[i] = dx0[i] + dy1[i];
668  }
669  return u;
670  }
671 
672  template<typename Scalar>
673  Func<Scalar>* init_fn(Solution<Scalar>*fu, const int order)
674  {
675  // Sanity checks.
676  if(fu == NULL) throw Hermes::Exceptions::Exception("NULL MeshFunction in Func<Scalar>*::init_fn().");
677  if(fu->get_mesh() == NULL) throw Hermes::Exceptions::Exception("Uninitialized MeshFunction used.");
678 
679  SpaceType space_type = fu->get_space_type();
680  SolutionType sln_type = fu->get_type();
681 
682  int nc = fu->get_num_components();
683  Quad2D* quad = fu->get_quad_2d();
684 #ifdef H2D_USE_SECOND_DERIVATIVES
685  if(space_type == HERMES_H1_SPACE && sln_type != HERMES_EXACT)
686  fu->set_quad_order(order, H2D_FN_ALL);
687  else
688  fu->set_quad_order(order);
689 #else
690  fu->set_quad_order(order);
691 #endif
692  double3* pt = quad->get_points(order, fu->get_active_element()->get_mode());
693  int np = quad->get_num_points(order, fu->get_active_element()->get_mode());
694  Func<Scalar>* u = new Func<Scalar>(np, nc);
695 
696  if(u->nc == 1)
697  {
698  u->val = new Scalar[np];
699  u->dx = new Scalar[np];
700  u->dy = new Scalar[np];
701 
702 #ifdef H2D_USE_SECOND_DERIVATIVES
703  if(space_type == HERMES_H1_SPACE && sln_type != HERMES_EXACT)
704  u->laplace = new Scalar[np];
705 #endif
706 
707  memcpy(u->val, fu->get_fn_values(), np * sizeof(Scalar));
708  memcpy(u->dx, fu->get_dx_values(), np * sizeof(Scalar));
709  memcpy(u->dy, fu->get_dy_values(), np * sizeof(Scalar));
710 
711 #ifdef H2D_USE_SECOND_DERIVATIVES
712  if(space_type == HERMES_H1_SPACE && sln_type != HERMES_EXACT)
713  {
714  Scalar *dxx = fu->get_dxx_values();
715  Scalar *dyy = fu->get_dyy_values();
716  for (int i = 0; i < np; i++)
717  u->laplace[i] = dxx[i] + dyy[i];
718  }
719 #endif
720  }
721  else if(u->nc == 2)
722  {
723  u->val0 = new Scalar[np];
724  u->val1 = new Scalar[np];
725  u->curl = new Scalar[np];
726  u->div = new Scalar[np];
727 
728  memcpy(u->val0, fu->get_fn_values(0), np * sizeof(Scalar));
729  memcpy(u->val1, fu->get_fn_values(1), np * sizeof(Scalar));
730 
731  Scalar *dx1 = fu->get_dx_values(1);
732  Scalar *dy0 = fu->get_dy_values(0);
733  for (int i = 0; i < np; i++)
734  u->curl[i] = dx1[i] - dy0[i];
735 
736  Scalar *dx0 = fu->get_dx_values(0);
737  Scalar *dy1 = fu->get_dy_values(1);
738  for (int i = 0; i < np; i++)
739  u->div[i] = dx0[i] + dy1[i];
740  }
741  return u;
742  }
743 
744  template Func<double>* init_fn(MeshFunction<double>*fu, const int order);
745  template Func<std::complex<double> >* init_fn(MeshFunction<std::complex<double> >*fu, const int order);
746 
747  template HERMES_API Func<double>* init_fn(Solution<double>*fu, const int order);
748  template HERMES_API Func<std::complex<double> >* init_fn(Solution<std::complex<double> >*fu, const int order);
749 
750  template class HERMES_API Func<Hermes::Ord>;
751  template class HERMES_API Func<double>;
752  template class HERMES_API Func<std::complex<double> >;
753  template class HERMES_API DiscontinuousFunc<Hermes::Ord>;
754  template class HERMES_API DiscontinuousFunc<double>;
755  template class HERMES_API DiscontinuousFunc<std::complex<double> >;
756  template class HERMES_API Geom<Hermes::Ord>;
757  template class HERMES_API Geom<double>;
758  template class HERMES_API InterfaceGeom<Hermes::Ord>;
759  template class HERMES_API InterfaceGeom<double>;
760  }
761 }