Hermes2D  2.0
solution.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 "exact_solution.h"
17 #include "forms.h"
18 
19 #include "solution_h2d_xml.h"
20 #include "ogprojection.h"
21 #include "api2d.h"
22 
23 #include <iostream>
24 #include <algorithm>
25 
26 namespace Hermes
27 {
28  namespace Hermes2D
29  {
30  static double3* cheb_tab_tri[11];
31  static double3* cheb_tab_quad[11];
32  static int cheb_np_tri[11];
33  static int cheb_np_quad[11];
34 
35  static double3** cheb_tab[2] = { cheb_tab_tri, cheb_tab_quad };
36  static int* cheb_np[2] = { cheb_np_tri, cheb_np_quad };
37 
38  static class Quad2DCheb : public Quad2D
39  {
40  public:
41 
42  Quad2DCheb()
43  {
44  max_order[0] = max_order[1] = 10;
45  num_tables[0] = num_tables[1] = 11;
46  tables = cheb_tab;
47  np = cheb_np;
48 
49  tables[0][0] = tables[1][0] = NULL;
50  np[0][0] = np[1][0] = 0;
51 
52  int i, j, k, n, m;
53  double3* pt;
54  for (int mode_i = 0; mode_i <= 1; mode_i++)
55  {
56  for (k = 0; k <= 10; k++)
57  {
58  np[mode_i][k] = n = mode_i ? sqr(k + 1) : (k + 1)*(k + 2)/2;
59  tables[mode_i][k] = pt = new double3[n];
60 
61  for (i = k, m = 0; i >= 0; i--)
62  for (j = k; j >= (mode_i ? 0 : k-i); j--, m++) {
63  pt[m][0] = k ? cos(j * M_PI / k) : 1.0;
64  pt[m][1] = k ? cos(i * M_PI / k) : 1.0;
65  pt[m][2] = 1.0;
66  }
67  }
68  }
69  };
70 
71  ~Quad2DCheb()
72  {
73  for (int mode_i = 0; mode_i <= 1; mode_i++)
74  for (int k = 1; k <= 10; k++)
75  delete [] tables[mode_i][k];
76  }
77 
78  virtual void dummy_fn() {}
79  } g_quad_2d_cheb;
80 
81  template<typename Scalar>
82  void Solution<Scalar>::set_static_verbose_output(bool verbose)
83  {
84  Solution<Scalar>::static_verbose_output = verbose;
85  }
86 
87  template<typename Scalar>
88  bool Solution<Scalar>::static_verbose_output = false;
89 
90  template<>
91  void Solution<double>::init()
92  {
93  memset(tables, 0, sizeof(tables));
94  memset(elems, 0, sizeof(elems));
95  memset(oldest, 0, sizeof(oldest));
96  transform = true;
97  sln_type = HERMES_UNDEF;
98  this->num_components = 0;
99  e_last = NULL;
100 
101  for(int i = 0; i < 4; i++)
102  for(int j = 0; j < 4; j++)
103  tables[i][j] = new std::map<uint64_t, LightArray<struct Function<double>::Node*>*>;
104 
105  mono_coeffs = NULL;
106  elem_coeffs[0] = elem_coeffs[1] = NULL;
107  elem_orders = NULL;
108  dxdy_buffer = NULL;
109  num_coeffs = num_elems = 0;
110  num_dofs = -1;
111 
112  this->set_quad_2d(&g_quad_2d_std);
113  }
114 
115  template<>
116  void Solution<std::complex<double> >::init()
117  {
118  memset(tables, 0, sizeof(tables));
119  memset(elems, 0, sizeof(elems));
120  memset(oldest, 0, sizeof(oldest));
121  transform = true;
122  sln_type = HERMES_UNDEF;
123  this->num_components = 0;
124  e_last = NULL;
125 
126  for(int i = 0; i < 4; i++)
127  for(int j = 0; j < 4; j++)
128  tables[i][j] = new std::map<uint64_t, LightArray<struct Function<std::complex<double> >::Node*>*>;
129 
130  mono_coeffs = NULL;
131  elem_coeffs[0] = elem_coeffs[1] = NULL;
132  elem_orders = NULL;
133  dxdy_buffer = NULL;
134  num_coeffs = num_elems = 0;
135  num_dofs = -1;
136 
137  this->set_quad_2d(&g_quad_2d_std);
138  }
139 
140  template<>
141  Solution<double>::Solution()
142  : MeshFunction<double>()
143  {
144  space_type = HERMES_INVALID_SPACE;
145  this->init();
146  }
147 
148  template<>
149  Solution<std::complex<double> >::Solution()
150  : MeshFunction<std::complex<double> >()
151  {
152  space_type = HERMES_INVALID_SPACE;
153  this->init();
154  }
155 
156  template<>
157  Solution<double>::Solution(const Mesh *mesh) : MeshFunction<double>(mesh)
158  {
159  space_type = HERMES_INVALID_SPACE;
160  this->init();
161  this->mesh = mesh;
162  }
163 
164  template<>
165  Solution<std::complex<double> >::Solution(const Mesh *mesh) : MeshFunction<std::complex<double> >(mesh)
166  {
167  space_type = HERMES_INVALID_SPACE;
168  this->init();
169  this->mesh = mesh;
170  }
171 
172  template<>
173  Solution<double>::Solution(Space<double>* s, Vector<double>* coeff_vec) : MeshFunction<double>(s->get_mesh())
174  {
175  space_type = s->get_type();
176  this->init();
177  this->mesh = s->get_mesh();
178  Solution<double>::vector_to_solution(coeff_vec, s, this);
179  }
180 
181  template<>
182  Solution<std::complex<double> >::Solution(Space<std::complex<double> >* s, Vector<std::complex<double> >* coeff_vec) : MeshFunction<std::complex<double> >(s->get_mesh())
183  {
184  space_type = s->get_type();
185  this->init();
186  this->mesh = s->get_mesh();
187  Solution<std::complex<double> >::vector_to_solution(coeff_vec, s, this);
188  }
189 
190  template<>
191  Solution<double>::Solution(Space<double>* s, double* coeff_vec) : MeshFunction<double>(s->get_mesh())
192  {
193  space_type = s->get_type();
194  this->init();
195  this->mesh = s->get_mesh();
196  Solution<double>::vector_to_solution(coeff_vec, s, this);
197  }
198 
199  template<>
200  Solution<std::complex<double> >::Solution(Space<std::complex<double> >* s, std::complex<double> * coeff_vec) : MeshFunction<std::complex<double> >(s->get_mesh())
201  {
202  space_type = s->get_type();
203  this->init();
204  this->mesh = s->get_mesh();
205  Solution<std::complex<double> >::vector_to_solution(coeff_vec, s, this);
206  }
207 
208  template<typename Scalar>
209  void Solution<Scalar>::assign(Solution<Scalar>* sln)
210  {
211  if(sln->sln_type == HERMES_UNDEF) throw Hermes::Exceptions::Exception("Solution being assigned is uninitialized.");
212  if(sln->sln_type != HERMES_SLN) { copy(sln); return; }
213 
214  free();
215 
216  this->mesh = sln->mesh;
217  // Solution vector and space setting.
218  space_type = sln->get_space_type();
219 
220  mono_coeffs = sln->mono_coeffs; sln->mono_coeffs = NULL;
221  elem_coeffs[0] = sln->elem_coeffs[0]; sln->elem_coeffs[0] = NULL;
222  elem_coeffs[1] = sln->elem_coeffs[1]; sln->elem_coeffs[1] = NULL;
223  elem_orders = sln->elem_orders; sln->elem_orders = NULL;
224  dxdy_buffer = sln->dxdy_buffer; sln->dxdy_buffer = NULL;
225  num_coeffs = sln->num_coeffs; sln->num_coeffs = 0;
226  num_elems = sln->num_elems; sln->num_elems = 0;
227 
228  sln_type = sln->sln_type;
229  this->num_components = sln->num_components;
230 
231  memset(sln->tables, 0, sizeof(sln->tables));
232  }
233 
234  template<typename Scalar>
235  void Solution<Scalar>::copy(const Solution<Scalar>* sln)
236  {
237  if(sln->sln_type == HERMES_UNDEF) throw Hermes::Exceptions::Exception("Solution being copied is uninitialized.");
238  free();
239 
240  this->mesh = sln->mesh;
241 
242  sln_type = sln->sln_type;
243  space_type = sln->get_space_type();
244  this->num_components = sln->num_components;
245  num_dofs = sln->num_dofs;
246 
247  if(sln->sln_type == HERMES_SLN) // standard solution: copy coefficient arrays
248  {
249  num_coeffs = sln->num_coeffs;
250  num_elems = sln->num_elems;
251 
252  mono_coeffs = new Scalar[num_coeffs];
253  memcpy(mono_coeffs, sln->mono_coeffs, sizeof(Scalar) * num_coeffs);
254 
255  for (int l = 0; l < this->num_components; l++)
256  {
257  elem_coeffs[l] = new int[num_elems];
258  memcpy(elem_coeffs[l], sln->elem_coeffs[l], sizeof(int) * num_elems);
259  }
260 
261  elem_orders = new int[num_elems];
262  memcpy(elem_orders, sln->elem_orders, sizeof(int) * num_elems);
263 
264  init_dxdy_buffer();
265  }
266  else // Const, exact handled differently.
267  throw Hermes::Exceptions::Exception("Undefined or exact solutions cannot be copied into an instance of Solution already coming from computation.");
268 
269  this->element = NULL;
270  }
271 
272  template<typename Scalar>
273  MeshFunction<Scalar>* Solution<Scalar>::clone() const
274  {
275  Solution<Scalar>* sln = new Solution<Scalar>();
276  sln->copy(this);
277  return sln;
278  }
279 
280  template<typename Scalar>
281  void Solution<Scalar>::free_tables()
282  {
283  for (int i = 0; i < 4; i++)
284  for (int j = 0; j < 4; j++)
285  if(tables[i][j] != NULL)
286  {
287  for(typename std::map<uint64_t, LightArray<struct Function<Scalar>::Node*>*>::iterator it = tables[i][j]->begin(); it != tables[i][j]->end(); it++)
288  {
289  for(unsigned int l = 0; l < it->second->get_size(); l++)
290  if(it->second->present(l))
291  ::free(it->second->get(l));
292  delete it->second;
293  }
294  tables[i][j]->clear();
295  delete tables[i][j];
296  tables[i][j] = NULL;
297  elems[i][j] = NULL;
298  }
299  }
300 
301  template<>
302  void Solution<double>::free()
303  {
304  if(mono_coeffs != NULL) { delete [] mono_coeffs; mono_coeffs = NULL; }
305  if(elem_orders != NULL) { delete [] elem_orders; elem_orders = NULL; }
306  if(dxdy_buffer != NULL) { delete [] dxdy_buffer; dxdy_buffer = NULL; }
307 
308  for (int i = 0; i < this->num_components; i++)
309  if(elem_coeffs[i] != NULL)
310  { delete [] elem_coeffs[i]; elem_coeffs[i] = NULL; }
311 
312  e_last = NULL;
313 
314  free_tables();
315  }
316 
317  template<>
318  void Solution<std::complex<double> >::free()
319  {
320  if(mono_coeffs != NULL) { delete [] mono_coeffs; mono_coeffs = NULL; }
321  if(elem_orders != NULL) { delete [] elem_orders; elem_orders = NULL; }
322  if(dxdy_buffer != NULL) { delete [] dxdy_buffer; dxdy_buffer = NULL; }
323 
324  for (int i = 0; i < this->num_components; i++)
325  if(elem_coeffs[i] != NULL)
326  { delete [] elem_coeffs[i]; elem_coeffs[i] = NULL; }
327 
328  e_last = NULL;
329 
330  free_tables();
331  }
332 
333  template<>
334  Solution<double>::~Solution()
335  {
336  free();
337  space_type = HERMES_INVALID_SPACE;
338  }
339 
340  template<>
341  Solution<std::complex<double> >::~Solution()
342  {
343  free();
344  space_type = HERMES_INVALID_SPACE;
345  }
346 
347  static struct mono_lu_init
348  {
349  public:
350 
351  // this is a set of LU-decomposed matrices shared by all Solutions
352  double** mat[2][11];
353  int* perm[2][11];
354 
355  mono_lu_init()
356  {
357  memset(mat, 0, sizeof(mat));
358  }
359 
360  ~mono_lu_init()
361  {
362  for (int m = 0; m <= 1; m++)
363  for (int i = 0; i <= 10; i++)
364  if(mat[m][i] != NULL) {
365  delete [] mat[m][i];
366  delete [] perm[m][i];
367  }
368  }
369  }
370  mono_lu;
371 
372  template<typename Scalar>
373  double** Solution<Scalar>::calc_mono_matrix(int o, int*& perm)
374  {
375  int i, j, k, l, m, row;
376  double x, y, xn, yn;
377  int n = this->mode ? sqr(o + 1) : (o + 1)*(o + 2)/2;
378 
379  // loop through all chebyshev points
380  double** mat = new_matrix<double>(n, n);
381  for (k = o, row = 0; k >= 0; k--)
382  {
383  y = o ? cos(k * M_PI / o) : 1.0;
384  for (l = o; l >= (this->mode ? 0 : o-k); l--, row++)
385  {
386  x = o ? cos(l * M_PI / o) : 1.0;
387 
388  // each row of the matrix contains all the monomials x^i*y^j
389  for (i = 0, yn = 1.0, m = n-1; i <= o; i++, yn *= y)
390  for (j = (this->mode ? 0 : i), xn = 1.0; j <= o; j++, xn *= x, m--)
391  mat[row][m] = xn * yn;
392  }
393  }
394 
395  double d;
396  perm = new int[n];
397  ludcmp(mat, n, perm, &d);
398  return mat;
399  }
400 
401  template<typename Scalar>
402  void Solution<Scalar>::set_coeff_vector(const Space<Scalar>* space, const Vector<Scalar>* vec,
403  bool add_dir_lift, int start_index)
404  {
405  // Sanity check.
406  if(space == NULL) throw Exceptions::NullException(1);
407  if(vec == NULL) throw Exceptions::NullException(2);
408 
409  space_type = space->get_type();
410  Scalar* coeffs = new Scalar[vec->length()];
411  vec->extract(coeffs);
412  this->set_coeff_vector(space, coeffs, add_dir_lift, start_index);
413  delete [] coeffs;
414  }
415 
416  template<typename Scalar>
417  void Solution<Scalar>::set_coeff_vector(const Space<Scalar>* space, const Scalar* coeffs,
418  bool add_dir_lift, int start_index)
419  {
420  // Sanity check.
421  if(space == NULL) throw Exceptions::NullException(1);
422 
423  // Initialize precalc shapeset using the space's shapeset.
424  Shapeset *shapeset = space->shapeset;
425  if(space->shapeset == NULL) throw Exceptions::Exception("Space->shapeset == NULL in Solution<Scalar>::set_coeff_vector().");
426  PrecalcShapeset *pss = new PrecalcShapeset(shapeset);
427  if(pss == NULL) throw Exceptions::Exception("PrecalcShapeset could not be allocated in Solution<Scalar>::set_coeff_vector().");
428  set_coeff_vector(space, pss, coeffs, add_dir_lift, start_index);
429  delete pss;
430  }
431 
432  template<typename Scalar>
433  void Solution<Scalar>::set_coeff_vector(const Space<Scalar>* space, PrecalcShapeset* pss,
434  const Scalar* coeff_vec, bool add_dir_lift, int start_index)
435  {
436  int o;
437  if(Solution<Scalar>::static_verbose_output)
438  Hermes::Mixins::Loggable::Static::info("Solution: set_coeff_vector called.");
439  // Sanity checks.
440  if(space == NULL) throw Exceptions::NullException(1);
441  if(space->get_mesh() == NULL) throw Exceptions::Exception("Mesh == NULL in Solution<Scalar>::set_coeff_vector().");
442  if(pss == NULL) throw Exceptions::NullException(2);
443  if(coeff_vec == NULL) throw Exceptions::NullException(3);
444  if(coeff_vec == NULL) throw Exceptions::Exception("Coefficient vector == NULL in Solution<Scalar>::set_coeff_vector().");
445  if(!space->is_up_to_date())
446  throw Exceptions::Exception("Provided 'space' is not up to date.");
447  if(space->shapeset != pss->shapeset)
448  throw Exceptions::Exception("Provided 'space' and 'pss' must have the same shapesets.");
449 
450  if(Solution<Scalar>::static_verbose_output)
451  Hermes::Mixins::Loggable::Static::info("Solution: set_coeff_vector - solution being freed.");
452 
453  free();
454 
455  this->space_type = space->get_type();
456 
457  this->num_components = pss->get_num_components();
458  sln_type = HERMES_SLN;
459 
460  // Copy the mesh.
461  this->mesh = space->get_mesh();
462 
463  // Allocate the coefficient arrays.
464  num_elems = this->mesh->get_max_element_id();
465  if(elem_orders != NULL)
466  delete [] elem_orders;
467  elem_orders = new int[num_elems];
468  memset(elem_orders, 0, sizeof(int) * num_elems);
469  for (int l = 0; l < this->num_components; l++)
470  {
471  if(elem_coeffs[l] != NULL)
472  delete [] elem_coeffs[l];
473  elem_coeffs[l] = new int[num_elems];
474  memset(elem_coeffs[l], 0, sizeof(int) * num_elems);
475  }
476 
477  // Obtain element orders, allocate mono_coeffs.
478  Element* e;
479  num_coeffs = 0;
480  for_all_active_elements(e, this->mesh)
481  {
482  this->mode = e->get_mode();
483  o = space->get_element_order(e->id);
484  o = std::max(H2D_GET_H_ORDER(o), H2D_GET_V_ORDER(o));
485  for (unsigned int k = 0; k < e->get_nvert(); k++)
486  {
487  int eo = space->get_edge_order(e, k);
488  if(eo > o) o = eo;
489  }
490 
491  // Hcurl and Hdiv: actual order of functions is one higher than element order
492  if((space->shapeset)->get_num_components() == 2) o++;
493 
494  num_coeffs += this->mode ? sqr(o + 1) : (o + 1)*(o + 2)/2;
495  elem_orders[e->id] = o;
496  }
497  num_coeffs *= this->num_components;
498  if(mono_coeffs != NULL)
499  delete [] mono_coeffs;
500  mono_coeffs = new Scalar[num_coeffs];
501 
502  // Express the solution on elements as a linear combination of monomials.
503  Quad2D* quad = &g_quad_2d_cheb;
504  pss->set_quad_2d(quad);
505  Scalar* mono = mono_coeffs;
506  for_all_active_elements(e, this->mesh)
507  {
508  this->mode = e->get_mode();
509  o = elem_orders[e->id];
510  int np = quad->get_num_points(o, e->get_mode());
511 
512  AsmList<Scalar> al;
513  space->get_element_assembly_list(e, &al);
514  pss->set_active_element(e);
515 
516  for (int l = 0; l < this->num_components; l++)
517  {
518  // Obtain solution values for the current element.
519  Scalar* val = mono;
520  elem_coeffs[l][e->id] = (int) (mono - mono_coeffs);
521  memset(val, 0, sizeof(Scalar)*np);
522  for (unsigned int k = 0; k < al.cnt; k++)
523  {
524  pss->set_active_shape(al.idx[k]);
525  pss->set_quad_order(o, H2D_FN_VAL);
526  int dof = al.dof[k];
527  double dir_lift_coeff = add_dir_lift ? 1.0 : 0.0;
528  // By subtracting space->first_dof we make sure that it does not matter where the
529  // enumeration of dofs in the space starts. This ca be either zero or there can be some
530  // offset. By adding start_index we move to the desired section of coeff_vec.
531  Scalar coef = al.coef[k] * (dof >= 0 ? coeff_vec[dof - space->first_dof + start_index] : dir_lift_coeff);
532  double* shape = pss->get_fn_values(l);
533  for (int i = 0; i < np; i++)
534  val[i] += shape[i] * coef;
535  }
536  mono += np;
537 
538  // solve for the monomial coefficients
539  if(mono_lu.mat[this->mode][o] == NULL)
540  mono_lu.mat[this->mode][o] = calc_mono_matrix(o, mono_lu.perm[this->mode][o]);
541  lubksb(mono_lu.mat[this->mode][o], np, mono_lu.perm[this->mode][o], val);
542  }
543  }
544 
545  if(this->mesh == NULL) throw Hermes::Exceptions::Exception("mesh == NULL.\n");
546  init_dxdy_buffer();
547  this->element = NULL;
548  if(Solution<Scalar>::static_verbose_output)
549  Hermes::Mixins::Loggable::Static::info("Solution: set_coeff_vector - done.");
550  }
551 
552  template<typename Scalar>
553  void Solution<Scalar>::vector_to_solutions(const Scalar* solution_vector,
554  Hermes::vector<const Space<Scalar>*> spaces, Hermes::vector<Solution<Scalar>*> solutions,
555  Hermes::vector<bool> add_dir_lift, Hermes::vector<int> start_indices)
556  {
557  if(solution_vector == NULL) throw Exceptions::NullException(1);
558  if(spaces.size() != solutions.size()) throw Exceptions::LengthException(2, 3, spaces.size(), solutions.size());
559 
560  // If start indices are not given, calculate them using the dimension of each space.
561  Hermes::vector<int> start_indices_new;
562  if(start_indices.empty())
563  {
564  int counter = 0;
565  for (int i=0; i < spaces.size(); i++)
566  {
567  start_indices_new.push_back(counter);
568  counter += spaces[i]->get_num_dofs();
569  }
570  }
571  else
572  {
573  if(start_indices.size() != spaces.size()) throw Hermes::Exceptions::Exception("Mismatched start indices in vector_to_solutions().");
574  for (int i=0; i < spaces.size(); i++)
575  {
576  start_indices_new.push_back(start_indices[i]);
577  }
578  }
579 
580  for(unsigned int i = 0; i < solutions.size(); i++)
581  {
582  if(Solution<Scalar>::static_verbose_output)
583  {
584  if(i == 0)
585  Hermes::Mixins::Loggable::Static::info("Solution: %d-th solution", i);
586  if(i == 1)
587  Hermes::Mixins::Loggable::Static::info("Solution: %d-st solution", i);
588  if(i == 2)
589  Hermes::Mixins::Loggable::Static::info("Solution: %d-nd solution", i);
590  if(i == 3)
591  Hermes::Mixins::Loggable::Static::info("Solution: %d-rd solution", i);
592  if(i > 3)
593  Hermes::Mixins::Loggable::Static::info("Solution: %d-th solution", i);
594  }
595  if(add_dir_lift == Hermes::vector<bool>())
596  Solution<Scalar>::vector_to_solution(solution_vector, spaces[i], solutions[i], true, start_indices_new[i]);
597  else
598  Solution<Scalar>::vector_to_solution(solution_vector, spaces[i], solutions[i], add_dir_lift[i], start_indices_new[i]);
599  }
600 
601  return;
602  }
603 
604  template<typename Scalar>
605  void Solution<Scalar>::vector_to_solution(const Scalar* solution_vector, const Space<Scalar>* space,
606  Solution<Scalar>* solution, bool add_dir_lift, int start_index)
607  {
608  // Sanity checks.
609  if(solution_vector == NULL) throw Exceptions::NullException(1);
610  if(space == NULL) throw Exceptions::NullException(2);
611  if(solution == NULL) throw Exceptions::NullException(3);
612 
613  solution->set_coeff_vector(space, solution_vector, add_dir_lift, start_index);
614  }
615 
616  template<typename Scalar>
617  void Solution<Scalar>::vector_to_solutions(const Vector<Scalar>* solution_vector, Hermes::vector<const Space<Scalar>*> spaces,
618  Hermes::vector<Solution<Scalar>*> solutions, Hermes::vector<bool> add_dir_lift, Hermes::vector<int> start_indices)
619  {
620  if(solution_vector == NULL) throw Exceptions::NullException(1);
621  if(spaces.size() != solutions.size()) throw Exceptions::LengthException(2, 3, spaces.size(), solutions.size());
622 
623  // If start indices are not given, calculate them using the dimension of each space.
624  Hermes::vector<int> start_indices_new;
625  if(start_indices.empty())
626  {
627  int counter = 0;
628  for (int i=0; i < spaces.size(); i++)
629  {
630  start_indices_new.push_back(counter);
631  counter += spaces[i]->get_num_dofs();
632  }
633  }
634  else
635  {
636  if(start_indices.size() != spaces.size()) throw Hermes::Exceptions::Exception("Mismatched start indices in vector_to_solutions().");
637  for (int i=0; i < spaces.size(); i++)
638  {
639  start_indices_new.push_back(start_indices[i]);
640  }
641  }
642 
643  for(unsigned int i = 0; i < solutions.size(); i++)
644  {
645  if(Solution<Scalar>::static_verbose_output)
646  {
647  if(i == 0)
648  Hermes::Mixins::Loggable::Static::info("Solution: %d-th solution", i);
649  if(i == 1)
650  Hermes::Mixins::Loggable::Static::info("Solution: %d-st solution", i);
651  if(i == 2)
652  Hermes::Mixins::Loggable::Static::info("Solution: %d-nd solution", i);
653  if(i == 3)
654  Hermes::Mixins::Loggable::Static::info("Solution: %d-rd solution", i);
655  if(i > 3)
656  Hermes::Mixins::Loggable::Static::info("Solution: %d-th solution", i);
657  }
658  if(add_dir_lift == Hermes::vector<bool>())
659  Solution<Scalar>::vector_to_solution(solution_vector, spaces[i], solutions[i], true, start_indices_new[i]);
660  else
661  Solution<Scalar>::vector_to_solution(solution_vector, spaces[i], solutions[i], add_dir_lift[i], start_indices_new[i]);
662  }
663 
664  return;
665  }
666 
667  template<typename Scalar>
668  void Solution<Scalar>::vector_to_solutions_common_dir_lift(const Vector<Scalar>* solution_vector, Hermes::vector<const Space<Scalar>*> spaces,
669  Hermes::vector<Solution<Scalar>*> solutions, bool add_dir_lift)
670  {
671  if(solution_vector == NULL) throw Exceptions::NullException(1);
672  if(spaces.size() != solutions.size()) throw Exceptions::LengthException(2, 3, spaces.size(), solutions.size());
673 
674  // If start indices are not given, calculate them using the dimension of each space.
675  Hermes::vector<int> start_indices_new;
676  int counter = 0;
677  for (int i=0; i < spaces.size(); i++)
678  {
679  start_indices_new.push_back(counter);
680  counter += spaces[i]->get_num_dofs();
681  }
682 
683  for(unsigned int i = 0; i < solutions.size(); i++)
684  Solution<Scalar>::vector_to_solution(solution_vector, spaces[i], solutions[i], add_dir_lift, start_indices_new[i]);
685 
686  return;
687  }
688 
689  template<typename Scalar>
690  void Solution<Scalar>::vector_to_solutions_common_dir_lift(const Scalar* solution_vector, Hermes::vector<const Space<Scalar>*> spaces,
691  Hermes::vector<Solution<Scalar>*> solutions, bool add_dir_lift)
692  {
693  if(solution_vector == NULL) throw Exceptions::NullException(1);
694  if(spaces.size() != solutions.size()) throw Exceptions::LengthException(2, 3, spaces.size(), solutions.size());
695 
696  // If start indices are not given, calculate them using the dimension of each space.
697  Hermes::vector<int> start_indices_new;
698  int counter = 0;
699  for (int i=0; i < spaces.size(); i++)
700  {
701  start_indices_new.push_back(counter);
702  counter += spaces[i]->get_num_dofs();
703  }
704 
705  for(unsigned int i = 0; i < solutions.size(); i++)
706  Solution<Scalar>::vector_to_solution(solution_vector, spaces[i], solutions[i], add_dir_lift, start_indices_new[i]);
707 
708  return;
709  }
710 
711  template<typename Scalar>
712  void Solution<Scalar>::vector_to_solution(const Vector<Scalar>* solution_vector, const Space<Scalar>* space,
713  Solution<Scalar>* solution, bool add_dir_lift, int start_index)
714  {
715  // Sanity checks.
716  if(solution_vector == NULL) throw Exceptions::NullException(1);
717  if(space == NULL) throw Exceptions::NullException(2);
718  if(solution == NULL) throw Exceptions::NullException(3);
719 
720  solution->set_coeff_vector(space, solution_vector, add_dir_lift, start_index);
721  }
722 
723  template<typename Scalar>
724  void Solution<Scalar>::vector_to_solutions(const Scalar* solution_vector, Hermes::vector<const Space<Scalar>*> spaces,
725  Hermes::vector<Solution<Scalar>*> solutions, Hermes::vector<PrecalcShapeset *> pss,
726  Hermes::vector<bool> add_dir_lift, Hermes::vector<int> start_indices)
727  {
728  if(solution_vector==NULL) throw Exceptions::NullException(1);
729  if(spaces.size() != solutions.size()) throw Exceptions::LengthException(2, 3, spaces.size(), solutions.size());
730 
731  // If start indices are not given, calculate them using the dimension of each space.
732  Hermes::vector<int> start_indices_new;
733  if(start_indices.empty())
734  {
735  int counter = 0;
736  for (int i=0; i < spaces.size(); i++)
737  {
738  start_indices_new.push_back(counter);
739  counter += spaces[i]->get_num_dofs();
740  }
741  }
742  else
743  {
744  if(start_indices.size() != spaces.size()) throw Hermes::Exceptions::Exception("Mismatched start indices in vector_to_solutions().");
745  for (int i=0; i < spaces.size(); i++)
746  {
747  start_indices_new.push_back(start_indices[i]);
748  }
749  }
750 
751  for(unsigned int i = 0; i < solutions.size(); i++)
752  {
753  if(Solution<Scalar>::static_verbose_output)
754  {
755  if(i == 0)
756  Hermes::Mixins::Loggable::Static::info("Solution: %d-th solution", i);
757  if(i == 1)
758  Hermes::Mixins::Loggable::Static::info("Solution: %d-st solution", i);
759  if(i == 2)
760  Hermes::Mixins::Loggable::Static::info("Solution: %d-nd solution", i);
761  if(i == 3)
762  Hermes::Mixins::Loggable::Static::info("Solution: %d-rd solution", i);
763  if(i > 3)
764  Hermes::Mixins::Loggable::Static::info("Solution: %d-th solution", i);
765  }
766  if(add_dir_lift == Hermes::vector<bool>())
767  Solution<Scalar>::vector_to_solution(solution_vector, spaces[i], solutions[i], true, start_indices_new[i]);
768  else
769  Solution<Scalar>::vector_to_solution(solution_vector, spaces[i], solutions[i], add_dir_lift[i], start_indices_new[i]);
770  }
771 
772  return;
773  }
774 
775  template<typename Scalar>
776  void Solution<Scalar>::vector_to_solution(const Scalar* solution_vector, const Space<Scalar>* space, Solution<Scalar>* solution,
777  PrecalcShapeset* pss, bool add_dir_lift, int start_index)
778  {
779  if(solution_vector == NULL) throw Exceptions::NullException(1);
780  if(space == NULL) throw Exceptions::NullException(2);
781  if(solution == NULL) throw Exceptions::NullException(3);
782  if(pss == NULL) throw Exceptions::NullException(4);
783 
784  solution->set_coeff_vector(space, pss, solution_vector, add_dir_lift, start_index);
785  }
786 
787  template<typename Scalar>
788  void Solution<Scalar>::set_dirichlet_lift(const Space<Scalar>* space, PrecalcShapeset* pss)
789  {
790  space_type = space->get_type();
791  int ndof = space->get_num_dofs();
792  Scalar *temp = new Scalar[ndof];
793  memset(temp, 0, sizeof(Scalar)*ndof);
794  bool add_dir_lift = true;
795  int start_index = 0;
796  this->set_coeff_vector(space, pss, temp, add_dir_lift, start_index);
797  delete [] temp;
798  }
799 
800  template<typename Scalar>
801  void Solution<Scalar>::enable_transform(bool enable)
802  {
803  if(transform != enable) free_tables();
804  transform = enable;
805  }
806 
807  template<typename Scalar>
808  void Solution<Scalar>::multiply(Scalar coef)
809  {
810  if(sln_type == HERMES_SLN)
811  {
812  for (int i = 0; i < num_coeffs; i++)
813  mono_coeffs[i] *= coef;
814  }
815  else if(sln_type == HERMES_EXACT)
816  dynamic_cast<ExactSolution<Scalar>*>(this)->exact_multiplicator *= coef;
817  else
818  throw Hermes::Exceptions::Exception("Uninitialized solution.");
819  }
820 
821  template<typename Scalar>
822  static void make_dx_coeffs(int mode, int o, Scalar* mono, Scalar* result)
823  {
824  int i, j, k;
825  for (i = 0; i <= o; i++) {
826  *result++= 0.0;
827  k = mode ? o : i;
828  for (j = 0; j < k; j++)
829  *result++= (Scalar) (k-j) * mono[j];
830  mono += k + 1;
831  }
832  }
833 
834  template<typename Scalar>
835  static void make_dy_coeffs(int mode, int o, Scalar* mono, Scalar* result)
836  {
837  int i, j;
838  if(mode) {
839  for (j = 0; j <= o; j++)
840  *result++= 0.0;
841  for (i = 0; i < o; i++)
842  for (j = 0; j <= o; j++)
843  *result++= (Scalar) (o-i) * (*mono++);
844  }
845  else {
846  for (i = 0; i <= o; i++) {
847  *result++= 0.0;
848  for (j = 0; j < i; j++)
849  *result++= (Scalar) (o + 1-i) * (*mono++);
850  }
851  }
852  }
853 
854  template<typename Scalar>
855  void Solution<Scalar>::init_dxdy_buffer()
856  {
857  if(dxdy_buffer != NULL)
858  {
859  delete [] dxdy_buffer;
860  dxdy_buffer = NULL;
861  }
862  dxdy_buffer = new Scalar[this->num_components * 5 * 121];
863  }
864 
865  template<typename Scalar>
866  void Solution<Scalar>::set_active_element(Element* e)
867  {
868  if(e == this->element)
869  return; // FIXME
870  if(!e->active) throw Hermes::Exceptions::Exception("Cannot select inactive element. Wrong mesh?");
871  MeshFunction<Scalar>::set_active_element(e);
872 
873  // try finding an existing table for e
874  for (cur_elem = 0; cur_elem < 4; cur_elem++)
875  if(elems[this->cur_quad][cur_elem] == e)
876  break;
877 
878  // if not found, free the oldest one and use its slot
879  if(cur_elem >= 4)
880  {
881  if(tables[this->cur_quad][oldest[this->cur_quad]] != NULL)
882  {
883  for(typename std::map<uint64_t, LightArray<struct Function<Scalar>::Node*>*>::iterator it = tables[this->cur_quad][oldest[this->cur_quad]]->begin(); it != tables[this->cur_quad][oldest[this->cur_quad]]->end(); it++)
884  {
885  for(unsigned int l = 0; l < it->second->get_size(); l++)
886  if(it->second->present(l))
887  ::free(it->second->get(l));
888  delete it->second;
889  }
890  delete tables[this->cur_quad][oldest[this->cur_quad]];
891  tables[this->cur_quad][oldest[this->cur_quad]] = NULL;
892  elems[this->cur_quad][oldest[this->cur_quad]] = NULL;
893  }
894 
895  tables[this->cur_quad][oldest[this->cur_quad]] = new std::map<uint64_t, LightArray<struct Function<Scalar>::Node*>*>;
896 
897  cur_elem = oldest[this->cur_quad];
898  if(++oldest[this->cur_quad] >= 4)
899  oldest[this->cur_quad] = 0;
900 
901  elems[this->cur_quad][cur_elem] = e;
902  }
903 
904  if(sln_type == HERMES_SLN)
905  {
906  int o = this->order = elem_orders[this->element->id];
907  int n = this->mode ? sqr(o + 1) : (o + 1)*(o + 2)/2;
908 
909  for (int i = 0, m = 0; i < this->num_components; i++)
910  {
911  Scalar* mono = mono_coeffs + elem_coeffs[i][e->id];
912  dxdy_coeffs[i][0] = mono;
913 
914  make_dx_coeffs(this->mode, o, mono, dxdy_coeffs[i][1] = dxdy_buffer + m); m += n;
915  make_dy_coeffs(this->mode, o, mono, dxdy_coeffs[i][2] = dxdy_buffer + m); m += n;
916  make_dx_coeffs(this->mode, o, dxdy_coeffs[i][1], dxdy_coeffs[i][3] = dxdy_buffer + m); m += n;
917  make_dy_coeffs(this->mode, o, dxdy_coeffs[i][2], dxdy_coeffs[i][4] = dxdy_buffer + m); m += n;
918  make_dx_coeffs(this->mode, o, dxdy_coeffs[i][2], dxdy_coeffs[i][5] = dxdy_buffer + m); m += n;
919  }
920  }
921  else if(sln_type == HERMES_EXACT)
922  {
923  this->order = Hermes::Hermes2D::g_max_quad;
924  }
925  else
926  throw Hermes::Exceptions::Exception("Uninitialized solution.");
927 
928  this->sub_tables = tables[this->cur_quad][cur_elem];
929 
930  this->update_nodes_ptr();
931  }
932 
933  template<typename Scalar>
934  static inline void set_vec_num(int n, Scalar* y, Scalar num)
935  {
936  for (int i = 0; i < n; i++)
937  y[i] = num;
938  }
939 
940  template<typename Scalar>
941  static inline void vec_x_vec_p_num(int n, Scalar* y, Scalar* x, Scalar num)
942  {
943  for (int i = 0; i < n; i++)
944  y[i] = y[i]*x[i] + num;
945  }
946 
947  template<typename Scalar>
948  static inline void vec_x_vec_p_vec(int n, Scalar* y, Scalar* x, Scalar* z)
949  {
950  for (int i = 0; i < n; i++)
951  y[i] = y[i]*x[i] + z[i];
952  }
953 
954  static const int H2D_GRAD = H2D_FN_DX_0 | H2D_FN_DY_0;
955  static const int H2D_SECOND = H2D_FN_DXX_0 | H2D_FN_DXY_0 | H2D_FN_DYY_0;
956  static const int H2D_CURL = H2D_FN_DX | H2D_FN_DY;
957 
958  template<typename Scalar>
959  void Solution<Scalar>::transform_values(int order, struct Function<Scalar>::Node* node, int newmask, int oldmask, int np)
960  {
961  double2x2 *mat, *m;
962  double3x2 *mat2, *mm;
963  int i, mstep = 0;
964 
965  // H1 space
966  if(space_type == HERMES_H1_SPACE)
967  {
968 #ifdef H2D_USE_SECOND_DERIVATIVES
969  if(((newmask & H2D_SECOND) == H2D_SECOND && (oldmask & H2D_SECOND) != H2D_SECOND))
970  {
971  this->update_refmap();
972  mat = this->refmap->get_inv_ref_map(order);
973  mat2 = this->refmap->get_second_ref_map(order);
974  for (i = 0, m = mat, mm = mat2; i < np; i++, m++, mm++)
975  {
976  Scalar vx = node->values[0][1][i];
977  Scalar vy = node->values[0][2][i];
978  Scalar vxx = node->values[0][3][i];
979  Scalar vyy = node->values[0][4][i];
980  Scalar vxy = node->values[0][5][i];
981 
982  node->values[0][3][i] = sqr((*m)[0][0])*vxx + 2*(*m)[0][1]*(*m)[0][0]*vxy + sqr((*m)[0][1])*vyy + (*mm)[0][0]*vx + (*mm)[0][1]*vy; // dxx
983  node->values[0][4][i] = sqr((*m)[1][0])*vxx + 2*(*m)[1][1]*(*m)[1][0]*vxy + sqr((*m)[1][1])*vyy + (*mm)[2][0]*vx + (*mm)[2][1]*vy; // dyy
984  node->values[0][5][i] = (*m)[0][0]*(*m)[1][0]*vxx + ((*m)[0][0]*(*m)[1][1] + (*m)[1][0]*(*m)[0][1])*vxy + (*m)[0][1]*(*m)[1][1]*vyy + (*mm)[1][0]*vx + (*mm)[1][1]*vy; //dxy
985  }
986  }
987 #endif
988  if((newmask & H2D_GRAD) == H2D_GRAD && (oldmask & H2D_GRAD) != H2D_GRAD)
989  {
990  this->update_refmap();
991  mat = this->refmap->get_const_inv_ref_map();
992  if(!this->refmap->is_jacobian_const()) { mat = this->refmap->get_inv_ref_map(order); mstep = 1; }
993 
994  for (i = 0, m = mat; i < np; i++, m += mstep)
995  {
996  Scalar vx = node->values[0][1][i];
997  Scalar vy = node->values[0][2][i];
998  node->values[0][1][i] = (*m)[0][0]*vx + (*m)[0][1]*vy;
999  node->values[0][2][i] = (*m)[1][0]*vx + (*m)[1][1]*vy;
1000  }
1001  }
1002  }
1003 
1004  // Hcurl space
1005  else if(space_type == HERMES_HCURL_SPACE)
1006  {
1007  bool trans_val = false, trans_curl = false;
1008  if((newmask & H2D_FN_VAL) == H2D_FN_VAL && (oldmask & H2D_FN_VAL) != H2D_FN_VAL) trans_val = true;
1009  if((newmask & H2D_CURL) == H2D_CURL && (oldmask & H2D_CURL) != H2D_CURL) trans_curl = true;
1010 
1011  if(trans_val || trans_curl)
1012  {
1013  this->update_refmap();
1014  mat = this->refmap->get_const_inv_ref_map();
1015  if(!this->refmap->is_jacobian_const()) { mat = this->refmap->get_inv_ref_map(order); mstep = 1; }
1016 
1017  for (i = 0, m = mat; i < np; i++, m += mstep)
1018  {
1019  if(trans_val)
1020  {
1021  Scalar vx = node->values[0][0][i];
1022  Scalar vy = node->values[1][0][i];
1023  node->values[0][0][i] = (*m)[0][0]*vx + (*m)[0][1]*vy;
1024  node->values[1][0][i] = (*m)[1][0]*vx + (*m)[1][1]*vy;
1025  }
1026  if(trans_curl)
1027  {
1028  Scalar e0x = node->values[0][1][i], e0y = node->values[0][2][i];
1029  Scalar e1x = node->values[1][1][i], e1y = node->values[1][2][i];
1030  node->values[1][1][i] = (*m)[0][0]*((*m)[1][0]*e0x + (*m)[1][1]*e1x) + (*m)[0][1]*((*m)[1][0]*e0y + (*m)[1][1]*e1y);
1031  node->values[0][2][i] = (*m)[1][0]*((*m)[0][0]*e0x + (*m)[0][1]*e1x) + (*m)[1][1]*((*m)[0][0]*e0y + (*m)[0][1]*e1y);
1032  }
1033  }
1034  }
1035  }
1036 
1037  // Hdiv space
1038  else if(space_type == HERMES_HDIV_SPACE)
1039  {
1040  if((newmask & H2D_FN_VAL) == H2D_FN_VAL && (oldmask & H2D_FN_VAL) != H2D_FN_VAL)
1041  {
1042  this->update_refmap();
1043  mat = this->refmap->get_const_inv_ref_map();
1044  if(!this->refmap->is_jacobian_const()) { mat = this->refmap->get_inv_ref_map(order); mstep = 1; }
1045 
1046  for (i = 0, m = mat; i < np; i++, m += mstep)
1047  {
1048  Scalar vx = node->values[0][0][i];
1049  Scalar vy = node->values[1][0][i];
1050  node->values[0][0][i] = (*m)[1][1]*vx - (*m)[1][0]*vy;
1051  node->values[1][0][i] = - (*m)[0][1]*vx + (*m)[0][0]*vy;
1052  }
1053  }
1054  }
1055  }
1056 
1057  template<typename Scalar>
1058  void Solution<Scalar>::precalculate(int order, int mask)
1059  {
1060  int i, j, k, l;
1061  struct Function<Scalar>::Node* node = NULL;
1062  Quad2D* quad = this->quads[this->cur_quad];
1063  int np = quad->get_num_points(order, this->mode);
1064 
1065  if(sln_type == HERMES_SLN)
1066  {
1067  // if we are required to transform vectors, we must precalculate both their components
1068  const int H2D_GRAD = H2D_FN_DX_0 | H2D_FN_DY_0;
1069  const int H2D_SECOND = H2D_FN_DXX_0 | H2D_FN_DXY_0 | H2D_FN_DYY_0;
1070  const int H2D_CURL = H2D_FN_DX | H2D_FN_DY; // sic
1071  if(transform)
1072  {
1073  if(this->num_components == 1) // H1 space or L2 space
1074  {
1075  if((mask & H2D_FN_DX_0) || (mask & H2D_FN_DY_0)) mask |= H2D_GRAD;
1076  if((mask & H2D_FN_DXX_0) || (mask & H2D_FN_DXY_0) || (mask & H2D_FN_DYY_0)) mask |= H2D_SECOND;
1077  }
1078  else if(space_type == HERMES_HCURL_SPACE) // Hcurl space
1079  { if((mask & H2D_FN_VAL_0) || (mask & H2D_FN_VAL_1)) mask |= H2D_FN_VAL;
1080  if((mask & H2D_FN_DX_1) || (mask & H2D_FN_DY_0)) mask |= H2D_CURL; }
1081  else // Hdiv space
1082  { if((mask & H2D_FN_VAL_0) || (mask & H2D_FN_VAL_1)) mask |= H2D_FN_VAL; }
1083  }
1084 
1085  int oldmask = (this->cur_node != NULL) ? this->cur_node->mask : 0;
1086  int newmask = mask | oldmask;
1087  node = this->new_node(newmask, np);
1088 
1089  // transform integration points by the current matrix
1090  Scalar* x = new Scalar[np];
1091  Scalar* y = new Scalar[np];
1092  Scalar* tx = new Scalar[np];
1093  double3* pt = quad->get_points(order, this->element->get_mode());
1094  for (i = 0; i < np; i++)
1095  {
1096  x[i] = pt[i][0] * this->ctm->m[0] + this->ctm->t[0];
1097  y[i] = pt[i][1] * this->ctm->m[1] + this->ctm->t[1];
1098  }
1099 
1100  // obtain the solution values, this is the core of the whole module
1101  int o = elem_orders[this->element->id];
1102  for (l = 0; l < this->num_components; l++)
1103  {
1104  for (k = 0; k < 6; k++)
1105  {
1106  if(newmask & this->idx2mask[k][l])
1107  {
1108  Scalar* result = node->values[l][k];
1109  if(oldmask & this->idx2mask[k][l])
1110  {
1111  // copy the old table if we have it already
1112  memcpy(result, this->cur_node->values[l][k], np * sizeof(Scalar));
1113  }
1114  else
1115  {
1116  // calculate the solution values using Horner's scheme
1117  Scalar* mono = dxdy_coeffs[l][k];
1118  for (i = 0; i <= o; i++)
1119  {
1120  set_vec_num(np, tx, *mono++);
1121  for (j = 1; j <= (this->mode ? o : i); j++)
1122  vec_x_vec_p_num(np, tx, x, *mono++);
1123 
1124  if(!i)
1125  memcpy(result, tx, sizeof(Scalar)*np);
1126  else
1127  vec_x_vec_p_vec(np, result, y, tx);
1128  }
1129  }
1130  }
1131  }
1132  }
1133 
1134  delete [] x;
1135  delete [] y;
1136  delete [] tx;
1137 
1138  // transform gradient or vector solution, if required
1139  if(transform)
1140  transform_values(order, node, newmask, oldmask, np);
1141  }
1142  else if(sln_type == HERMES_EXACT)
1143  {
1144  if(mask & ~H2D_FN_DEFAULT)
1145  throw Hermes::Exceptions::Exception("Cannot obtain second derivatives of an exact solution.");
1146  node = this->new_node(mask = H2D_FN_DEFAULT, np);
1147 
1148  this->update_refmap();
1149  double* x = this->refmap->get_phys_x(order);
1150  double* y = this->refmap->get_phys_y(order);
1151 
1152  // evaluate the exact solution
1153  if(this->num_components == 1)
1154  {
1155  // untransform values
1156  if(!transform)
1157  {
1158  double2x2 *mat, *m;
1159  int mstep = 0;
1160  mat = this->refmap->get_const_inv_ref_map();
1161  if(!this->refmap->is_jacobian_const()) { mat = this->refmap->get_inv_ref_map(order); mstep = 1; }
1162 
1163  for (i = 0, m = mat; i < np; i++, m += mstep)
1164  {
1165  double jac = (*m)[0][0] * (*m)[1][1] - (*m)[1][0] * (*m)[0][1];
1166  Scalar val, dx = 0.0, dy = 0.0;
1167  val = (static_cast<ExactSolutionScalar<Scalar>*>(this))->exact_function(x[i], y[i], dx, dy);
1168  node->values[0][0][i] = val * (static_cast<ExactSolutionVector<Scalar>*>(this))->exact_multiplicator;
1169  node->values[0][1][i] = ( (*m)[1][1]*dx - (*m)[0][1]*dy) / jac * (static_cast<ExactSolutionVector<Scalar>*>(this))->exact_multiplicator;
1170  node->values[0][2][i] = (- (*m)[1][0]*dx + (*m)[0][0]*dy) / jac * (static_cast<ExactSolutionVector<Scalar>*>(this))->exact_multiplicator;
1171  }
1172  }
1173  else
1174  {
1175  for (i = 0; i < np; i++)
1176  {
1177  Scalar val, dx = 0.0, dy = 0.0;
1178  val = (static_cast<ExactSolutionScalar<Scalar>*>(this))->exact_function(x[i], y[i], dx, dy);
1179  node->values[0][0][i] = val * (static_cast<ExactSolutionVector<Scalar>*>(this))->exact_multiplicator;
1180  node->values[0][1][i] = dx * (static_cast<ExactSolutionVector<Scalar>*>(this))->exact_multiplicator;
1181  node->values[0][2][i] = dy * (static_cast<ExactSolutionVector<Scalar>*>(this))->exact_multiplicator;
1182  }
1183  }
1184  }
1185  else
1186  {
1187  for (i = 0; i < np; i++)
1188  {
1189  Scalar2<Scalar> dx (0.0, 0.0 ), dy ( 0.0, 0.0 );
1190  Scalar2<Scalar> val = (static_cast<ExactSolutionVector<Scalar>*>(this))->exact_function(x[i], y[i], dx, dy);
1191  for (j = 0; j < 2; j++)
1192  {
1193  node->values[j][0][i] = val[j] * (static_cast<ExactSolutionVector<Scalar>*>(this))->exact_multiplicator;
1194  node->values[j][1][i] = dx[j] * (static_cast<ExactSolutionVector<Scalar>*>(this))->exact_multiplicator;
1195  node->values[j][2][i] = dy[j] * (static_cast<ExactSolutionVector<Scalar>*>(this))->exact_multiplicator;
1196  }
1197  }
1198  }
1199  }
1200  else
1201  {
1202  throw Hermes::Exceptions::Exception("Cannot obtain values -- uninitialized solution. The solution was either "
1203  "not calculated yet or you used the assignment operator which destroys "
1204  "the solution on its right-hand side.");
1205  }
1206 
1207  if(this->nodes->present(order))
1208  {
1209  assert(this->nodes->get(order) == this->cur_node);
1210  ::free(this->nodes->get(order));
1211  }
1212  this->nodes->add(node, order);
1213  this->cur_node = node;
1214  }
1215 
1216  template<>
1217  void Solution<double>::save(const char* filename) const
1218  {
1219  if(sln_type == HERMES_UNDEF)
1220  throw Exceptions::Exception("Cannot save -- uninitialized solution.");
1221 
1222  try
1223  {
1224  XMLSolution::solution xmlsolution(this->num_components, this->num_elems, this->num_coeffs, 0, 0);
1225  switch(this->get_space_type())
1226  {
1227  case HERMES_H1_SPACE:
1228  xmlsolution.space().set("h1");
1229  break;
1230  case HERMES_HCURL_SPACE:
1231  xmlsolution.space().set("hcurl");
1232  break;
1233  case HERMES_HDIV_SPACE:
1234  xmlsolution.space().set("hdiv");
1235  break;
1236  case HERMES_L2_SPACE:
1237  xmlsolution.space().set("l2");
1238  break;
1239  }
1240 
1241  for(unsigned int coeffs_i = 0; coeffs_i < this->num_coeffs; coeffs_i++)
1242  xmlsolution.mono_coeffs().push_back(XMLSolution::mono_coeffs(coeffs_i, mono_coeffs[coeffs_i]));
1243 
1244  for(unsigned int elems_i = 0; elems_i < this->num_elems; elems_i++)
1245  xmlsolution.elem_orders().push_back(XMLSolution::elem_orders(elems_i, elem_orders[elems_i]));
1246 
1247  for (unsigned int component_i = 0; component_i < this->num_components; component_i++)
1248  {
1249  xmlsolution.component().push_back(XMLSolution::component());
1250  xmlsolution.component().back().component_number() = component_i;
1251  for(unsigned int elems_i = 0; elems_i < this->num_elems; elems_i++)
1252  xmlsolution.component().back().elem_coeffs().push_back(XMLSolution::elem_coeffs(elems_i, elem_coeffs[component_i][elems_i]));
1253  }
1254 
1255  std::string solution_schema_location(Hermes2DApi.get_text_param_value(xmlSchemasDirPath));
1256  solution_schema_location.append("/solution_h2d_xml.xsd");
1257  ::xml_schema::namespace_info namespace_info_solution("XMLSolution", solution_schema_location);
1258 
1259  ::xml_schema::namespace_infomap namespace_info_map;
1260  namespace_info_map.insert(std::pair<std::basic_string<char>, xml_schema::namespace_info>("solution", namespace_info_solution));
1261 
1262  std::ofstream out(filename);
1263 
1264  ::xml_schema::flags parsing_flags = ::xml_schema::flags::dont_pretty_print;
1265 
1266  XMLSolution::solution_(out, xmlsolution, namespace_info_map, "UTF-8", parsing_flags);
1267  out.close();
1268  }
1269  catch (const xml_schema::exception& e)
1270  {
1271  throw Hermes::Exceptions::SolutionSaveFailureException(e.what());
1272  }
1273  return;
1274  }
1275 
1276  template<>
1277  void Solution<std::complex<double> >::save(const char* filename) const
1278  {
1279  if(sln_type == HERMES_UNDEF)
1280  throw Exceptions::Exception("Cannot save -- uninitialized solution.");
1281 
1282  try
1283  {
1284  XMLSolution::solution xmlsolution(this->num_components, this->num_elems, this->num_coeffs, 0, 1);
1285 
1286  switch(this->get_space_type())
1287  {
1288  case HERMES_H1_SPACE:
1289  xmlsolution.space().set("h1");
1290  break;
1291  case HERMES_HCURL_SPACE:
1292  xmlsolution.space().set("hcurl");
1293  break;
1294  case HERMES_HDIV_SPACE:
1295  xmlsolution.space().set("hdiv");
1296  break;
1297  case HERMES_L2_SPACE:
1298  xmlsolution.space().set("l2");
1299  break;
1300  }
1301 
1302  for(unsigned int coeffs_i = 0; coeffs_i < this->num_coeffs; coeffs_i++)
1303  {
1304  xmlsolution.mono_coeffs().push_back(XMLSolution::mono_coeffs(coeffs_i, mono_coeffs[coeffs_i].real()));
1305  xmlsolution.mono_coeffs().back().im() = mono_coeffs[coeffs_i].imag();
1306  }
1307 
1308  for(unsigned int elems_i = 0; elems_i < this->num_elems; elems_i++)
1309  xmlsolution.elem_orders().push_back(XMLSolution::elem_orders(elems_i, elem_orders[elems_i]));
1310 
1311  for (unsigned int component_i = 0; component_i < this->num_components; component_i++)
1312  {
1313  xmlsolution.component().push_back(XMLSolution::component());
1314  xmlsolution.component().back().component_number() = component_i;
1315  for(unsigned int elems_i = 0; elems_i < this->num_elems; elems_i++)
1316  xmlsolution.component().back().elem_coeffs().push_back(XMLSolution::elem_coeffs(elems_i, elem_coeffs[component_i][elems_i]));
1317  }
1318 
1319  std::string solution_schema_location(Hermes2DApi.get_text_param_value(xmlSchemasDirPath));
1320  solution_schema_location.append("/solution_h2d_xml.xsd");
1321  ::xml_schema::namespace_info namespace_info_solution("XMLSolution", solution_schema_location);
1322 
1323  ::xml_schema::namespace_infomap namespace_info_map;
1324  namespace_info_map.insert(std::pair<std::basic_string<char>, xml_schema::namespace_info>("solution", namespace_info_solution));
1325 
1326  std::ofstream out(filename);
1327 
1328  ::xml_schema::flags parsing_flags = ::xml_schema::flags::dont_pretty_print;
1329 
1330  XMLSolution::solution_(out, xmlsolution, namespace_info_map, "UTF-8", parsing_flags);
1331  out.close();
1332  }
1333  catch (const xml_schema::exception& e)
1334  {
1335  throw Hermes::Exceptions::SolutionSaveFailureException(e.what());
1336  }
1337  return;
1338  }
1339 
1340  template<>
1341  void Solution<double>::load(const char* filename, Space<double>* space)
1342  {
1343  free();
1344  this->mesh = space->get_mesh();
1345  this->space_type = space->get_type();
1346 
1347  try
1348  {
1349  ::xml_schema::flags parsing_flags = 0;
1350  if(!this->validate)
1351  parsing_flags = xml_schema::flags::dont_validate;
1352 
1353  std::auto_ptr<XMLSolution::solution> parsed_xml_solution(XMLSolution::solution_(filename, parsing_flags));
1354  sln_type = parsed_xml_solution->exact() == 0 ? HERMES_SLN : HERMES_EXACT;
1355 
1356  if(sln_type == HERMES_EXACT)
1357  {
1358  switch(parsed_xml_solution->ncmp())
1359  {
1360  case 1:
1361  if(parsed_xml_solution->exactC() == 0)
1362  {
1363  double* coeff_vec = new double[space->get_num_dofs()];
1364  ConstantSolution<double> sln(mesh, parsed_xml_solution->exactCXR().get());
1365  OGProjection<double> ogProj;
1366  ogProj.project_global(space, &sln, coeff_vec);
1367  this->set_coeff_vector(space, coeff_vec, true, 0);
1368  sln_type = HERMES_SLN;
1369  }
1370  else
1371  throw Hermes::Exceptions::SolutionLoadFailureException("Mismatched real - complex exact solutions.");
1372  break;
1373  case 2:
1374  if(parsed_xml_solution->exactC() == 0)
1375  {
1376  double* coeff_vec = new double[space->get_num_dofs()];
1377  ConstantSolutionVector<double> sln(mesh, parsed_xml_solution->exactCXR().get(), parsed_xml_solution->exactCYR().get());
1378  OGProjection<double> ogProj;
1379  ogProj.project_global(space, &sln, coeff_vec);
1380  this->set_coeff_vector(space, coeff_vec, true, 0);
1381  sln_type = HERMES_SLN;
1382  }
1383  else
1384  throw Hermes::Exceptions::SolutionLoadFailureException("Mismatched real - complex exact solutions.");
1385  break;
1386  }
1387  }
1388  else
1389  {
1390  if(!strcmp(parsed_xml_solution->space().get().c_str(),"h1"))
1391  if(this->space_type != HERMES_H1_SPACE)
1392  throw Exceptions::Exception("Space types not compliant in Solution::load().");
1393 
1394  if(!strcmp(parsed_xml_solution->space().get().c_str(),"l2"))
1395  if(this->space_type != HERMES_L2_SPACE)
1396  throw Exceptions::Exception("Space types not compliant in Solution::load().");
1397 
1398  if(!strcmp(parsed_xml_solution->space().get().c_str(),"hcurl"))
1399  if(this->space_type != HERMES_HCURL_SPACE)
1400  throw Exceptions::Exception("Space types not compliant in Solution::load().");
1401 
1402  if(!strcmp(parsed_xml_solution->space().get().c_str(),"hdiv"))
1403  if(this->space_type != HERMES_HDIV_SPACE)
1404  throw Exceptions::Exception("Space types not compliant in Solution::load().");
1405 
1406  this->num_coeffs = parsed_xml_solution->nc();
1407  this->num_elems = parsed_xml_solution->nel();
1408  this->num_components = parsed_xml_solution->ncmp();
1409 
1410  this->mono_coeffs = new double[num_coeffs];
1411  memset(this->mono_coeffs, 0, this->num_coeffs*sizeof(double));
1412 
1413  for(unsigned int component_i = 0; component_i < num_components; component_i++)
1414  elem_coeffs[component_i] = new int[num_elems];
1415 
1416  this->elem_orders = new int[num_elems];
1417 
1418  for (unsigned int coeffs_i = 0; coeffs_i < num_coeffs; coeffs_i++)
1419  this->mono_coeffs[parsed_xml_solution->mono_coeffs().at(coeffs_i).id()] = parsed_xml_solution->mono_coeffs().at(coeffs_i).re();
1420 
1421  for (unsigned int elems_i = 0; elems_i < num_elems; elems_i++)
1422  this->elem_orders[parsed_xml_solution->elem_orders().at(elems_i).id()] = parsed_xml_solution->elem_orders().at(elems_i).ord();
1423 
1424  for (unsigned int component_i = 0; component_i < this->num_components; component_i++)
1425  for (unsigned int elems_i = 0; elems_i < num_elems; elems_i++)
1426  this->elem_coeffs[component_i][parsed_xml_solution->component().at(component_i).elem_coeffs().at(elems_i).id()] = parsed_xml_solution->component().at(component_i).elem_coeffs().at(elems_i).c();
1427 
1428  }
1429  init_dxdy_buffer();
1430  }
1431  catch (const xml_schema::exception& e)
1432  {
1433  throw Hermes::Exceptions::SolutionLoadFailureException(e.what());
1434  }
1435  return;
1436  }
1437 
1438  template<>
1439  void Solution<std::complex<double> >::load(const char* filename, Space<std::complex<double> >* space)
1440  {
1441  free();
1442  sln_type = HERMES_SLN;
1443  this->mesh = space->get_mesh();
1444  this->space_type = space->get_type();
1445 
1446  try
1447  {
1448  ::xml_schema::flags parsing_flags = 0;
1449  if(!this->validate)
1450  parsing_flags = xml_schema::flags::dont_validate;
1451 
1452  std::auto_ptr<XMLSolution::solution> parsed_xml_solution(XMLSolution::solution_(filename, parsing_flags));
1453  sln_type = parsed_xml_solution->exact() == 0 ? HERMES_SLN : HERMES_EXACT;
1454 
1455  if(sln_type == HERMES_EXACT)
1456  {
1457  switch(parsed_xml_solution->ncmp())
1458  {
1459  case 1:
1460  if(parsed_xml_solution->exactC() == 1)
1461  {
1462  std::complex<double>* coeff_vec = new std::complex<double>[space->get_num_dofs()];
1463  ConstantSolution<std::complex<double> > sln(mesh, std::complex<double>(parsed_xml_solution->exactCXR().get(), parsed_xml_solution->exactCXC().get()));
1464  OGProjection<std::complex<double> > ogProj;
1465  ogProj.project_global(space, &sln, coeff_vec);
1466  this->set_coeff_vector(space, coeff_vec, true, 0);
1467  sln_type = HERMES_SLN;
1468  }
1469  else
1470  throw Hermes::Exceptions::SolutionLoadFailureException("Mismatched real - complex exact solutions.");
1471  break;
1472  case 2:
1473  if(parsed_xml_solution->exactC() == 1)
1474  {
1475  std::complex<double>* coeff_vec = new std::complex<double>[space->get_num_dofs()];
1476  ConstantSolutionVector<std::complex<double> > sln(mesh, std::complex<double>(parsed_xml_solution->exactCXR().get(), parsed_xml_solution->exactCXC().get()), std::complex<double>(parsed_xml_solution->exactCYR().get(), parsed_xml_solution->exactCYC().get()));
1477  OGProjection<std::complex<double> > ogProj;
1478  ogProj.project_global(space, &sln, coeff_vec);
1479  this->set_coeff_vector(space, coeff_vec, true, 0);
1480  sln_type = HERMES_SLN;
1481  }
1482  else
1483  throw Hermes::Exceptions::SolutionLoadFailureException("Mismatched real - complex exact solutions.");
1484  break;
1485  }
1486  }
1487  else
1488  {
1489  if(!strcmp(parsed_xml_solution->space().get().c_str(),"h1"))
1490  if(this->space_type != HERMES_H1_SPACE)
1491  throw Exceptions::Exception("Space types not compliant in Solution::load().");
1492 
1493  if(!strcmp(parsed_xml_solution->space().get().c_str(),"l2"))
1494  if(this->space_type != HERMES_L2_SPACE)
1495  throw Exceptions::Exception("Space types not compliant in Solution::load().");
1496 
1497  if(!strcmp(parsed_xml_solution->space().get().c_str(),"hcurl"))
1498  if(this->space_type != HERMES_HCURL_SPACE)
1499  throw Exceptions::Exception("Space types not compliant in Solution::load().");
1500 
1501  if(!strcmp(parsed_xml_solution->space().get().c_str(),"hdiv"))
1502  if(this->space_type != HERMES_HDIV_SPACE)
1503  throw Exceptions::Exception("Space types not compliant in Solution::load().");
1504 
1505  ::xml_schema::flags parsing_flags = 0;
1506  if(!this->validate)
1507  parsing_flags = xml_schema::flags::dont_validate;
1508 
1509  std::auto_ptr<XMLSolution::solution> parsed_xml_solution(XMLSolution::solution_(filename, parsing_flags));
1510 
1511  this->num_coeffs = parsed_xml_solution->nc();
1512  this->num_elems = parsed_xml_solution->nel();
1513  this->num_components = parsed_xml_solution->ncmp();
1514 
1515  this->mono_coeffs = new std::complex<double>[num_coeffs];
1516  memset(this->mono_coeffs, 0, this->num_coeffs*sizeof(std::complex<double>));
1517 
1518  for(unsigned int component_i = 0; component_i < num_components; component_i++)
1519  elem_coeffs[component_i] = new int[num_elems];
1520 
1521  this->elem_orders = new int[num_elems];
1522 
1523  for (unsigned int coeffs_i = 0; coeffs_i < num_coeffs; coeffs_i++)
1524  this->mono_coeffs[parsed_xml_solution->mono_coeffs().at(coeffs_i).id()] = std::complex<double>(parsed_xml_solution->mono_coeffs().at(coeffs_i).re(), parsed_xml_solution->mono_coeffs().at(coeffs_i).im().get());
1525 
1526  for (unsigned int elems_i = 0; elems_i < num_elems; elems_i++)
1527  this->elem_orders[parsed_xml_solution->elem_orders().at(elems_i).id()] = parsed_xml_solution->elem_orders().at(elems_i).ord();
1528 
1529  for (unsigned int component_i = 0; component_i < this->num_components; component_i++)
1530  for (unsigned int elems_i = 0; elems_i < num_elems; elems_i++)
1531  this->elem_coeffs[component_i][parsed_xml_solution->component().at(component_i).elem_coeffs().at(elems_i).id()] = parsed_xml_solution->component().at(component_i).elem_coeffs().at(elems_i).c();
1532  }
1533 
1534  init_dxdy_buffer();
1535  }
1536  catch (const xml_schema::exception& e)
1537  {
1538  throw Hermes::Exceptions::SolutionLoadFailureException(e.what());
1539  }
1540  return;
1541  }
1542 
1543  template<typename Scalar>
1544  Scalar Solution<Scalar>::get_ref_value(Element* e, double xi1, double xi2, int component, int item)
1545  {
1546  if(e==NULL)
1547  throw Exceptions::NullException(1);
1548 
1549  set_active_element(e);
1550 
1551  int o = elem_orders[e->id];
1552  Scalar* mono = dxdy_coeffs[component][item];
1553  Scalar result = 0.0;
1554  int k = 0;
1555  for (int i = 0; i <= o; i++)
1556  {
1557  Scalar row = mono[k++];
1558  for (int j = 0; j < (this->mode ? o : i); j++)
1559  row = row * xi1 + mono[k++];
1560  result = result * xi2 + row;
1561  }
1562  return result;
1563  }
1564 
1565  template<typename Scalar>
1566  Scalar Solution<Scalar>::get_ref_value_transformed(Element* e, double xi1, double xi2, int a, int b)
1567  {
1568  if(e==NULL) throw Exceptions::NullException(1);
1569  if(this->num_components == 1)
1570  {
1571  if(b == 0)
1572  return get_ref_value(e, xi1, xi2, a, b);
1573  if(b == 1 || b == 2)
1574  {
1575  double2x2 m;
1576  double xx, yy;
1577  this->refmap->inv_ref_map_at_point(xi1, xi2, xx, yy, m);
1578  Scalar dx = get_ref_value(e_last = e, xi1, xi2, a, 1);
1579  Scalar dy = get_ref_value(e, xi1, xi2, a, 2);
1580  if(b == 1) return m[0][0]*dx + m[0][1]*dy; // H2D_FN_DX
1581  if(b == 2) return m[1][0]*dx + m[1][1]*dy; // H2D_FN_DY
1582  }
1583  else
1584  {
1585 #ifdef H2D_USE_SECOND_DERIVATIVES
1586  double2x2 mat;
1587  double3x2 mat2;
1588  double xx, yy;
1589 
1590  this->refmap->inv_ref_map_at_point(xi1, xi2, xx, yy, mat);
1591  this->refmap->second_ref_map_at_point(xi1, xi2, xx, yy, mat2);
1592 
1593  Scalar vx = get_ref_value(e, xi1, xi2, a, 1);
1594  Scalar vy = get_ref_value(e, xi1, xi2, a, 2);
1595  Scalar vxx = get_ref_value(e, xi1, xi2, a, 3);
1596  Scalar vyy = get_ref_value(e, xi1, xi2, a, 4);
1597  Scalar vxy = get_ref_value(e, xi1, xi2, a, 5);
1598  if(b == 3)
1599  return sqr(mat[0][0])*vxx + 2*mat[0][1]*mat[0][0]*vxy + sqr(mat[0][1])*vyy + mat2[0][0]*vx + mat2[0][1]*vy; // dxx
1600  if(b == 4)
1601  return sqr(mat[1][0])*vxx + 2*mat[1][1]*mat[1][0]*vxy + sqr(mat[1][1])*vyy + mat2[2][0]*vx + mat2[2][1]*vy; // dyy
1602  if(b == 5)
1603  return mat[0][0]*mat[1][0]*vxx + (mat[0][0]*mat[1][1] + mat[1][0]*mat[0][1])*vxy + mat[0][1]*mat[1][1]*vyy + mat2[1][0]*vx + mat2[1][1]*vy; //dxy
1604 #endif
1605  }
1606  }
1607  else // vector solution
1608  {
1609  if(b == 0)
1610  {
1611  double2x2 m;
1612  double xx, yy;
1613  this->refmap->inv_ref_map_at_point(xi1, xi2, xx, yy, m);
1614  Scalar vx = get_ref_value(e, xi1, xi2, 0, 0);
1615  Scalar vy = get_ref_value(e, xi1, xi2, 1, 0);
1616  if(a == 0) return m[0][0]*vx + m[0][1]*vy; // H2D_FN_VAL_0
1617  if(a == 1) return m[1][0]*vx + m[1][1]*vy; // H2D_FN_VAL_1
1618  }
1619  else
1620  throw Hermes::Exceptions::Exception("Getting derivatives of the vector solution: Not implemented yet.");
1621  }
1622  throw Hermes::Exceptions::Exception("internal error: reached end of non-void function");
1623  return 0;
1624  }
1625 
1626  template<typename Scalar>
1627  Func<Scalar>* Solution<Scalar>::get_pt_value(double x, double y)
1628  {
1629  double xi1, xi2;
1630 
1631  Func<Scalar>* toReturn = new Func<Scalar>(1, this->num_components);
1632 
1633  if(sln_type == HERMES_EXACT)
1634  {
1635  if(this->num_components == 1)
1636  {
1637  Scalar val, dx = 0.0, dy = 0.0;
1638  val = (static_cast<ExactSolutionScalar<Scalar>*>(this))->exact_function(x, y, dx, dy);
1639  toReturn->val = new Scalar[1];
1640  toReturn->dx = new Scalar[1];
1641  toReturn->dy = new Scalar[1];
1642  toReturn->val[0] = val * (static_cast<ExactSolutionScalar<Scalar>*>(this))->exact_multiplicator;
1643  toReturn->dx[0] = dx * (static_cast<ExactSolutionScalar<Scalar>*>(this))->exact_multiplicator;
1644  toReturn->dy[0] = dy * (static_cast<ExactSolutionScalar<Scalar>*>(this))->exact_multiplicator;
1645  }
1646  else
1647  {
1648  Scalar2<Scalar> dx(0.0, 0.0), dy(0.0, 0.0);
1649  Scalar2<Scalar> val = (static_cast<ExactSolutionVector<Scalar>*>(this))->exact_function(x, y, dx, dy);
1650 
1651  toReturn->val0 = new Scalar[1];
1652  toReturn->dx0 = new Scalar[1];
1653  toReturn->dy0 = new Scalar[1];
1654  toReturn->val1 = new Scalar[1];
1655  toReturn->dx1 = new Scalar[1];
1656  toReturn->dy1 = new Scalar[1];
1657  toReturn->val0[0] = val[0] * (static_cast<ExactSolutionScalar<Scalar>*>(this))->exact_multiplicator;
1658  toReturn->val1[0] = val[1] * (static_cast<ExactSolutionScalar<Scalar>*>(this))->exact_multiplicator;
1659  toReturn->dx0[0] = dx[0] * (static_cast<ExactSolutionScalar<Scalar>*>(this))->exact_multiplicator;
1660  toReturn->dx1[0] = dx[1] * (static_cast<ExactSolutionScalar<Scalar>*>(this))->exact_multiplicator;
1661  toReturn->dy0[0] = dy[0] * (static_cast<ExactSolutionScalar<Scalar>*>(this))->exact_multiplicator;
1662  toReturn->dy1[0] = dy[1] * (static_cast<ExactSolutionScalar<Scalar>*>(this))->exact_multiplicator;
1663  }
1664 #ifdef H2D_USE_SECOND_DERIVATIVES
1665  this->warn("Cannot obtain second derivatives of an exact solution.");
1666 #endif
1667  }
1668  else if(sln_type == HERMES_UNDEF)
1669  {
1670  throw Hermes::Exceptions::Exception("Cannot obtain values -- uninitialized solution. The solution was either "
1671  "not calculated yet or you used the assignment operator which destroys "
1672  "the solution on its right-hand side.");
1673  }
1674  else // HERMES_SLN
1675  {
1676  Element* e = RefMap::element_on_physical_coordinates(this->mesh, x, y, &xi1, &xi2);
1677  if(e != NULL)
1678  {
1679  if(this->num_components == 1)
1680  {
1681  toReturn->val = new Scalar[1];
1682  toReturn->dx = new Scalar[1];
1683  toReturn->dy = new Scalar[1];
1684 
1685  toReturn->val[0] = get_ref_value(e, xi1, xi2, 0, 0);
1686 
1687  double2x2 m;
1688  double xx, yy;
1689  this->refmap->inv_ref_map_at_point(xi1, xi2, xx, yy, m);
1690  Scalar dx = get_ref_value(e, xi1, xi2, 0, 1);
1691  Scalar dy = get_ref_value(e, xi1, xi2, 0, 2);
1692  toReturn->dx[0] = m[0][0]*dx + m[0][1]*dy;
1693  toReturn->dy[0] = m[1][0]*dx + m[1][1]*dy;
1694 
1695 #ifdef H2D_USE_SECOND_DERIVATIVES
1696  toReturn->laplace = new Scalar[1];
1697  double2x2 mat;
1698  double3x2 mat2;
1699 
1700  this->refmap->inv_ref_map_at_point(xi1, xi2, xx, yy, mat);
1701  this->refmap->second_ref_map_at_point(xi1, xi2, xx, yy, mat2);
1702 
1703  Scalar vxx = get_ref_value(e, xi1, xi2, 0, 3);
1704  Scalar vyy = get_ref_value(e, xi1, xi2, 0, 4);
1705  Scalar vxy = get_ref_value(e, xi1, xi2, 0, 5);
1706  Scalar dxx = sqr(mat[0][0])*vxx + 2*mat[0][1]*mat[0][0]*vxy + sqr(mat[0][1])*vyy + mat2[0][0]*dx + mat2[0][1]*dy; // dxx
1707  Scalar dyy = sqr(mat[1][0])*vxx + 2*mat[1][1]*mat[1][0]*vxy + sqr(mat[1][1])*vyy + mat2[2][0]*dx + mat2[2][1]*dy; // dyy
1708  toReturn->laplace[0] = dxx + dyy;
1709 #endif
1710  }
1711  else // vector solution
1712  {
1713  toReturn->val0 = new Scalar[1];
1714  toReturn->val1 = new Scalar[1];
1715 
1716  double2x2 m;
1717  double xx, yy;
1718  this->refmap->inv_ref_map_at_point(xi1, xi2, xx, yy, m);
1719  Scalar vx = get_ref_value(e, xi1, xi2, 0, 0);
1720  Scalar vy = get_ref_value(e, xi1, xi2, 1, 0);
1721  toReturn->val0[0] = m[0][0]*vx + m[0][1]*vy;
1722  toReturn->val1[0] = m[1][0]*vx + m[1][1]*vy;
1723  Hermes::Mixins::Loggable::Static::warn("Derivatives of vector functions not implemented yet.");
1724  }
1725  return toReturn;
1726  }
1727 
1728  this->warn("Point (%g, %g) does not lie in any element.", x, y);
1729  return NULL;
1730  }
1731  }
1732 
1733  template class HERMES_API Solution<double>;
1734  template class HERMES_API Solution<std::complex<double> >;
1735  }
1736 }