Hermes2D  2.0
curved.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 "curved.h"
17 #include <algorithm>
18 #include "global.h"
19 #include "shapeset/shapeset_h1_all.h"
20 #include "shapeset/shapeset_common.h"
21 #include "shapeset/precalc.h"
22 #include "mesh.h"
23 #include "quad_all.h"
24 #include "matrix.h"
25 
26 using namespace Hermes::Algebra::DenseMatrixOperations;
27 namespace Hermes
28 {
29  namespace Hermes2D
30  {
31  double** CurvMap::edge_proj_matrix = NULL;
32  double** CurvMap::bubble_proj_matrix_tri = NULL;
33  double** CurvMap::bubble_proj_matrix_quad = NULL;
34 
35  double* CurvMap::edge_p = NULL;
36  double* CurvMap::bubble_tri_p = NULL;
37  double* CurvMap::bubble_quad_p = NULL;
38 
39  Quad1DStd CurvMap::quad1d;
40  Quad2DStd CurvMap::quad2d;
41 
42  Trf CurvMap::ctm;
43 
44  static double lambda_0(double x, double y)
45  {
46  return -0.5 * (x + y);
47  }
48  static double lambda_1(double x, double y)
49  {
50  return 0.5 * (x + 1);
51  }
52  static double lambda_2(double x, double y)
53  {
54  return 0.5 * (y + 1);
55  }
56 
57  bool CurvMap::warning_issued = false;
58 
59  double CurvMap::nurbs_basis_fn(int i, int k, double t, double* knot)
60  {
61  if(k == 0)
62  {
63  return (t >= knot[i] && t <= knot[i + 1] && knot[i] < knot[i + 1]) ? 1.0 : 0.0;
64  }
65  else
66  {
67  double N1 = nurbs_basis_fn(i, k-1, t, knot);
68  double N2 = nurbs_basis_fn(i + 1, k-1, t, knot);
69 
70  double result = 0.0;
71  if(knot[i + k] != knot[i])
72  {
73  result += ((t - knot[i]) / (knot[i + k] - knot[i])) * N1;
74  }
75  if(knot[i + k+1] != knot[i + 1])
76  {
77  result += ((knot[i + k+1] - t) / (knot[i + k+1] - knot[i + 1])) * N2;
78  }
79  return result;
80  }
81  }
82 
83  void CurvMap::nurbs_edge(Element* e, Nurbs* nurbs, int edge, double t, double& x,
84  double& y, double& n_x, double& n_y, double& t_x, double& t_y)
85  {
86  // Nurbs curves are parametrized from 0 to 1.
87  t = (t + 1.0) / 2.0;
88 
89  // Start point A, end point B.
90  double2 A, B;
91  A[0] = e->vn[edge]->x;
92  A[1] = e->vn[edge]->y;
93  B[0] = e->vn[e->next_vert(edge)]->x;
94  B[1] = e->vn[e->next_vert(edge)]->y;
95 
96  // Vector pointing from A to B.
97  double2 v;
98  v[0] = B[0] - A[0];
99  v[1] = B[1] - A[1];
100  double abs_v = sqrt(sqr(v[0]) + sqr(v[1]));
101 
102  // Straight line.
103  if(nurbs == NULL)
104  {
105  x = A[0] + t * v[0];
106  y = A[1] + t * v[1];
107  t_x = v[0] / abs_v;
108  t_y = v[1] / abs_v;
109  n_x = t_y;
110  n_y = -t_x;
111  }
112  else
113  {
114  // Circular arc.
115  if(nurbs->arc)
116  {
117  double3* cp = nurbs->pt;
118  x = y = 0.0;
119  double sum = 0.0; // sum of basis fns and weights
120 
121  for (int i = 0; i < nurbs->np; i++)
122  {
123  double basis = nurbs_basis_fn(i, nurbs->degree, t, nurbs->kv);
124  sum += cp[i][2] * basis;
125  double x_i = cp[i][0];
126  double y_i = cp[i][1];
127  double w_i = cp[i][2];
128  x += w_i * basis * x_i;
129  y += w_i * basis * y_i;
130  }
131 
132  x /= sum;
133  y /= sum;
134 
135  // Normal and tangential vectors.
136  // FIXME; This calculation is artificial and it assumes that
137  // the NURBS is a circular arc. This should be done in the
138  // same way for all NURBS.
139 
140  // End points, midpoint.
141  double2 M;
142  M[0] = (A[0] + B[0]) / 2.;
143  M[1] = (A[1] + B[1]) / 2.;
144  //printf("***** A = %g %g\n", A[0], A[1]);
145  //printf("***** B = %g %g\n", B[0], B[1]);
146  //printf("***** M = %g %g\n", M[0], M[1]);
147 
148  // Unit vector from M to center of circle S.
149  double2 w;
150  w[0] = -v[1] / abs_v;
151  w[1] = v[0] / abs_v;
152 
153  // Distance L between M and center of circle S
154  // can be calculated using a right-angle triangle
155  // whose one angle is alpha/2.
156  double alpha_rad = nurbs->angle * M_PI / 180.;
157  double L = 0.5 * abs_v / tan(0.5 * alpha_rad);
158  //printf("***** L = %g\n", L);
159 
160  // Center of circle.
161  double2 S;
162  S[0] = M[0] + w[0] * L;
163  S[1] = M[1] + w[1] * L;
164  //printf("***** S = %g %g\n", S[0], S[1]);
165 
166  // Calculation of radius and test.
167  double2 SA, SB;
168  SA[0] = A[0] - S[0];
169  SA[1] = A[1] - S[1];
170  SB[0] = B[0] - S[0];
171  SB[1] = B[1] - S[1];
172  double R = sqrt(sqr(SA[0]) + sqr(SA[1]));
173  double R2 = sqrt(sqr(SB[0]) + sqr(SB[1]));
174  if(std::abs(R - R2) > 1e-6)
175  throw Hermes::Exceptions::Exception("Internal error in nurbs_edge() - bad radius R.");
176 
177  // Normal vectors to circular arc at edge end points A, B.
178  double2 normal_A, normal_B;
179  normal_A[0] = SA[0] / R;
180  normal_A[1] = SA[1] / R;
181  normal_B[0] = SB[0] / R;
182  normal_B[1] = SB[1] / R;
183  //printf("***** normal_A = %g %g\n", normal_A[0], normal_A[1]);
184  //printf("***** normal_B = %g %g\n", normal_B[0], normal_B[1]);
185 
186  // Calculate rotational matrix that transforms AS_ref = (R, 0) -> SA
187  // and SB_ref -> SB.
188  double2 SB_ref;
189  SB_ref[0] = R * cos(alpha_rad);
190  SB_ref[1] = R * sin(alpha_rad);
191  // First we need to invert the matrix[(R 0)^T, SB_ref^T].
192  double m_11, m_12, m_21, m_22;
193  m_11 = R;
194  m_12 = SB_ref[0];
195  m_21 = 0;
196  m_22 = SB_ref[1];
197  double m_det = m_11 * m_22 - m_12 * m_21;
198  double inv_11, inv_12, inv_21, inv_22;
199  inv_11 = m_22 / m_det;
200  inv_12 = -m_12 / m_det;
201  inv_21 = -m_21 / m_det;
202  inv_22 = m_11 / m_det;
203  double s_11, s_12, s_21, s_22;
204  s_11 = SA[0];
205  s_12 = SB[0];
206  s_21 = SA[1];
207  s_22 = SB[1];
208  //Rotation matrix.
209  double r_11, r_12, r_21, r_22;
210  r_11 = s_11 * inv_11 + s_12 * inv_21;
211  r_12 = s_11 * inv_12 + s_12 * inv_22;
212  r_21 = s_21 * inv_11 + s_22 * inv_21;
213  r_22 = s_21 * inv_12 + s_22 * inv_22;
214  // The desired normal vector in reference coordinates.
215  double n_x_ref = cos(alpha_rad * t);
216  double n_y_ref = sin(alpha_rad * t);
217  // Rotate it.
218  n_x = r_11 * n_x_ref + r_12 * n_y_ref;
219  n_y = r_21 * n_x_ref + r_22 * n_y_ref;
220 
221  /*
222  // Calculate normal at point corresponding to the
223  // position of parameter 't' between 0 and 1.
224  n_x = normal_A[0] + t * (normal_B[0] - normal_A[0]);
225  n_y = normal_A[1] + t * (normal_B[1] - normal_A[1]);
226  double size_n = sqrt(sqr(n_x) + sqr(n_y));
227  n_x /= size_n;
228  n_y /= size_n;
229  */
230 
231  // Calculate tangential vectors.
232  t_x = -n_y;
233  t_y = n_x;
234 
235  // Correcting sign so that the normal points outside
236  // if the angle is negative.
237  if(nurbs->angle < 0)
238  {
239  n_x *= -1;
240  n_y *= -1;
241  t_x *= -1;
242  t_y *= -1;
243  }
244  }
245  // General NURBS.
246  // FIXME - calculation of normal and tangential vectors needs to be added.
247  else
248  {
249  double3* cp = nurbs->pt;
250  x = y = 0.0;
251  double sum = 0.0; // sum of basis fns and weights
252 
253  for (int i = 0; i < nurbs->np; i++)
254  {
255  double basis = nurbs_basis_fn(i, nurbs->degree, t, nurbs->kv);
256  sum += cp[i][2] * basis;
257  double x_i = cp[i][0];
258  double y_i = cp[i][1];
259  double w_i = cp[i][2];
260  x += w_i * basis * x_i;
261  y += w_i * basis * y_i;
262  }
263 
264  x /= sum;
265  y /= sum;
266 
267  if(!warning_issued)
268  {
269  printf("FIXME: IMPLEMENT CALCULATION OF n_x, n_y, t_x, t_y in nurbs_edge() !!!\n");
270  warning_issued = true;
271  }
272  n_x = 0;
273  n_y = 0;
274  t_x = 0;
275  t_y = 0;
276  }
277  }
278  }
279 
281  const double2 CurvMap::ref_vert[2][H2D_MAX_NUMBER_VERTICES] =
282  {
283  { { -1.0, -1.0 }, { 1.0, -1.0 }, { -1.0, 1.0 }, { 0.0, 0.0 } },
284  { { -1.0, -1.0 }, { 1.0, -1.0 }, { 1.0, 1.0 }, { -1.0, 1.0 } }
285  };
286 
287  // subtraction of straight edge and nurbs curve
288  void CurvMap::nurbs_edge_0(Element* e, Nurbs* nurbs, int edge, double t, double& x, double& y, double& n_x, double& n_y, double& t_x, double& t_y)
289  {
290  int va = edge;
291  int vb = e->next_vert(edge);
292  nurbs_edge(e, nurbs, edge, t, x, y, n_x, n_y, t_x, t_y);
293 
294  x -= 0.5 * ((1-t) * (e->vn[va]->x) + (1 + t) * (e->vn[vb]->x));
295  y -= 0.5 * ((1-t) * (e->vn[va]->y) + (1 + t) * (e->vn[vb]->y));
296 
297  double k = 4.0 / ((1-t) * (1 + t));
298  x *= k;
299  y *= k;
300  }
301 
302  // calculation of nonpolynomial reference mapping on curved element
303  void CurvMap::calc_ref_map_tri(Element* e, Nurbs** nurbs, double xi_1, double xi_2, double& x, double& y)
304  {
305  double fx, fy;
306  x = y = 0.0;
307 
308  for (unsigned int j = 0; j < e->get_nvert(); j++)
309  {
310  int va = j;
311  int vb = e->next_vert(j);
312  double l_a = 0;
313  double l_b = 0;
314  switch(va)
315  {
316  case 0:
317  l_a = lambda_0(xi_1, xi_2);
318  break;
319  case 1:
320  l_a = lambda_1(xi_1, xi_2);
321  break;
322  case 2:
323  l_a = lambda_2(xi_1, xi_2);
324  break;
325  }
326 
327  switch(vb)
328  {
329  case 0:
330  l_b = lambda_0(xi_1, xi_2);
331  break;
332  case 1:
333  l_b = lambda_1(xi_1, xi_2);
334  break;
335  case 2:
336  l_b = lambda_2(xi_1, xi_2);
337  break;
338  }
339 
340  // vertex part
341  x += e->vn[j]->x * l_a;
342  y += e->vn[j]->y * l_a;
343 
344  if(!(((ref_vert[0][va][0] == xi_1) && (ref_vert[0][va][1] == xi_2)) ||
345  ((ref_vert[0][vb][0] == xi_1) && (ref_vert[0][vb][1] == xi_2))))
346  {
347  // edge part
348  double t = l_b - l_a;
349  double n_x, n_y, t_x, t_y;
350  nurbs_edge_0(e, nurbs[j], j, t, fx, fy, n_x, n_y, t_x, t_y);
351  x += fx * l_a * l_b;
352  y += fy * l_a * l_b;
353  }
354  }
355  }
356 
357  void CurvMap::calc_ref_map_quad(Element* e, Nurbs** nurbs, double xi_1, double xi_2,
358  double& x, double& y)
359  {
360  double ex[H2D_MAX_NUMBER_EDGES], ey[H2D_MAX_NUMBER_EDGES];
361 
362  double n_x, n_y, t_x, t_y;
363  nurbs_edge(e, nurbs[0], 0, xi_1, ex[0], ey[0], n_x, n_y, t_x, t_y);
364  nurbs_edge(e, nurbs[1], 1, xi_2, ex[1], ey[1], n_x, n_y, t_x, t_y);
365  nurbs_edge(e, nurbs[2], 2, -xi_1, ex[2], ey[2], n_x, n_y, t_x, t_y);
366  nurbs_edge(e, nurbs[3], 3, -xi_2, ex[3], ey[3], n_x, n_y, t_x, t_y);
367 
368  x = (1-xi_2)/2.0 * ex[0] + (1 + xi_1)/2.0 * ex[1] +
369  (1 + xi_2)/2.0 * ex[2] + (1-xi_1)/2.0 * ex[3] -
370  (1-xi_1)*(1-xi_2)/4.0 * e->vn[0]->x - (1 + xi_1)*(1-xi_2)/4.0 * e->vn[1]->x -
371  (1 + xi_1)*(1 + xi_2)/4.0 * e->vn[2]->x - (1-xi_1)*(1 + xi_2)/4.0 * e->vn[3]->x;
372 
373  y = (1-xi_2)/2.0 * ey[0] + (1 + xi_1)/2.0 * ey[1] +
374  (1 + xi_2)/2.0 * ey[2] + (1-xi_1)/2.0 * ey[3] -
375  (1-xi_1)*(1-xi_2)/4.0 * e->vn[0]->y - (1 + xi_1)*(1-xi_2)/4.0 * e->vn[1]->y -
376  (1 + xi_1)*(1 + xi_2)/4.0 * e->vn[2]->y - (1-xi_1)*(1 + xi_2)/4.0 * e->vn[3]->y;
377  }
378 
379  void CurvMap::calc_ref_map(Element* e, Nurbs** nurbs, double xi_1, double xi_2, double2& f)
380  {
381  if(e->get_mode() == HERMES_MODE_QUAD)
382  calc_ref_map_quad(e, nurbs, xi_1, xi_2, f[0], f[1]);
383  else
384  calc_ref_map_tri(e, nurbs, xi_1, xi_2, f[0], f[1]);
385  }
386 
388 
389  // preparation of projection matrices, Cholesky factorization
390  void CurvMap::precalculate_cholesky_projection_matrix_edge(H1ShapesetJacobi* ref_map_shapeset, PrecalcShapeset* ref_map_pss)
391  {
392  int order = ref_map_shapeset->get_max_order();
393  int n = order - 1; // number of edge basis functions
394 
395  if(!edge_proj_matrix)
396  edge_proj_matrix = new_matrix<double>(n, n);
397 
398  // calculate projection matrix of maximum order
399  for (int i = 0; i < n; i++)
400  {
401  for (int j = i; j < n; j++)
402  {
403  int o = i + j + 4;
404  double2* pt = quad1d.get_points(o);
405  double val = 0.0;
406  for (int k = 0; k < quad1d.get_num_points(o); k++)
407  {
408  double fi = 0;
409  double fj = 0;
410  double x = pt[k][0];
411  switch(i + 2)
412  {
413  case 0:
414  fi = l0(x);
415  break;
416  case 1:
417  fi = l1(x);
418  break;
419  case 2:
420  fi = l2(x);
421  break;
422  case 3:
423  fi = l3(x);
424  break;
425  case 4:
426  fi = l4(x);
427  break;
428  case 5:
429  fi = l5(x);
430  break;
431  case 6:
432  fi = l6(x);
433  break;
434  case 7:
435  fi = l7(x);
436  break;
437  case 8:
438  fi = l8(x);
439  break;
440  case 9:
441  fi = l9(x);
442  break;
443  case 10:
444  fi = l10(x);
445  break;
446  case 11:
447  fi = l11(x);
448  break;
449  }
450  switch(j + 2)
451  {
452  case 0:
453  fj = l0(x);
454  break;
455  case 1:
456  fj = l1(x);
457  break;
458  case 2:
459  fj = l2(x);
460  break;
461  case 3:
462  fj = l3(x);
463  break;
464  case 4:
465  fj = l4(x);
466  break;
467  case 5:
468  fj = l5(x);
469  break;
470  case 6:
471  fj = l6(x);
472  break;
473  case 7:
474  fj = l7(x);
475  break;
476  case 8:
477  fj = l8(x);
478  break;
479  case 9:
480  fj = l9(x);
481  break;
482  case 10:
483  fj = l10(x);
484  break;
485  case 11:
486  fj = l11(x);
487  break;
488  }
489  val += pt[k][1] * (fi * fj);
490  }
491  edge_proj_matrix[i][j] = edge_proj_matrix[j][i] = val;
492  }
493  }
494 
495  // Cholesky factorization of the matrix
496  if(!edge_p)
497  edge_p = new double[n];
498  choldc(edge_proj_matrix, n, edge_p);
499  }
500 
501  // calculate the H1 seminorm products (\phi_i, \phi_j) for all 0 <= i, j < n, n is the number of bubble functions
502  double** CurvMap::calculate_bubble_projection_matrix(int nb, int* indices, H1ShapesetJacobi* ref_map_shapeset, PrecalcShapeset* ref_map_pss, ElementMode2D mode)
503  {
504  double** mat = new_matrix<double>(nb, nb);
505 
506  for (int i = 0; i < nb; i++)
507  {
508  for (int j = i; j < nb; j++)
509  {
510  int ii = indices[i], ij = indices[j];
511  int o = ref_map_shapeset->get_order(ii, mode) + ref_map_shapeset->get_order(ij, mode);
512  o = std::max(H2D_GET_V_ORDER(o), H2D_GET_H_ORDER(o));
513 
514  ref_map_pss->set_active_shape(ii);
515  ref_map_pss->set_quad_order(o, H2D_FN_VAL);
516  double* fni = ref_map_pss->get_fn_values();
517 
518  ref_map_pss->set_active_shape(ij);
519  ref_map_pss->set_quad_order(o, H2D_FN_VAL);
520  double* fnj = ref_map_pss->get_fn_values();
521 
522  double3* pt = quad2d.get_points(o, mode);
523  double val = 0.0;
524  for (int k = 0; k < quad2d.get_num_points(o, mode); k++)
525  val += pt[k][2] * (fni[k] * fnj[k]);
526 
527  mat[i][j] = mat[j][i] = val;
528  }
529  }
530 
531  return mat;
532  }
533 
534  void CurvMap::precalculate_cholesky_projection_matrices_bubble(H1ShapesetJacobi* ref_map_shapeset, PrecalcShapeset* ref_map_pss)
535  {
536  // *** triangles ***
537  int order = ref_map_shapeset->get_max_order();
538 
539  // calculate projection matrix of maximum order
540  if(ref_map_pss->get_active_element()->get_mode() == HERMES_MODE_TRIANGLE)
541  {
542  int nb = ref_map_shapeset->get_num_bubbles(order, HERMES_MODE_TRIANGLE);
543  int* indices = ref_map_shapeset->get_bubble_indices(order, HERMES_MODE_TRIANGLE);
544  bubble_proj_matrix_tri = calculate_bubble_projection_matrix(nb, indices, ref_map_shapeset, ref_map_pss, HERMES_MODE_TRIANGLE);
545 
546  // cholesky factorization of the matrix
547  bubble_tri_p = new double[nb];
548  choldc(bubble_proj_matrix_tri, nb, bubble_tri_p);
549  }
550 
551  // *** quads ***
552 
553  // calculate projection matrix of maximum order
554  if(ref_map_pss->get_active_element()->get_mode() == HERMES_MODE_QUAD)
555  {
556  order = H2D_MAKE_QUAD_ORDER(order, order);
557  int nb = ref_map_shapeset->get_num_bubbles(order, HERMES_MODE_QUAD);
558  int *indices = ref_map_shapeset->get_bubble_indices(order, HERMES_MODE_QUAD);
559 
560  bubble_proj_matrix_quad = calculate_bubble_projection_matrix(nb, indices, ref_map_shapeset, ref_map_pss, HERMES_MODE_QUAD);
561 
562  // cholesky factorization of the matrix
563  bubble_quad_p = new double[nb];
564  choldc(bubble_proj_matrix_quad, nb, bubble_quad_p);
565  }
566  }
567 
569 
570  // compute point (x, y) in reference element, edge vector (v1, v2)
571  void CurvMap::edge_coord(Element* e, int edge, double t, double2& x, double2& v)
572  {
573  int mode = e->get_mode();
574  double2 a, b;
575  a[0] = ctm.m[0] * ref_vert[mode][edge][0] + ctm.t[0];
576  a[1] = ctm.m[1] * ref_vert[mode][edge][1] + ctm.t[1];
577  b[0] = ctm.m[0] * ref_vert[mode][e->next_vert(edge)][0] + ctm.t[0];
578  b[1] = ctm.m[1] * ref_vert[mode][e->next_vert(edge)][1] + ctm.t[1];
579 
580  for (int i = 0; i < 2; i++)
581  {
582  v[i] = b[i] - a[i];
583  x[i] = a[i] + (t + 1.0)/2.0 * v[i];
584  }
585  double lenght = sqrt(v[0] * v[0] + v[1] * v[1]);
586  v[0] /= lenght; v[1] /= lenght;
587  }
588 
589  void CurvMap::calc_edge_projection(Element* e, int edge, Nurbs** nurbs, int order, double2* proj, H1ShapesetJacobi* ref_map_shapeset, PrecalcShapeset* ref_map_pss)
590  {
591  ref_map_pss->set_active_element(e);
592 
593  int i, j, k;
594  int mo1 = quad1d.get_max_order();
595  int np = quad1d.get_num_points(mo1);
596  int ne = order - 1;
597  int mode = e->get_mode();
598 
599  assert(np <= 15 && ne <= 10);
600  double2 fn[15];
601  double rhside[2][10];
602  memset(fn, 0, sizeof(double2) * np);
603  memset(rhside[0], 0, sizeof(double) * ne);
604  memset(rhside[1], 0, sizeof(double) * ne);
605 
606  double a_1, a_2, b_1, b_2;
607  a_1 = ctm.m[0] * ref_vert[mode][edge][0] + ctm.t[0];
608  a_2 = ctm.m[1] * ref_vert[mode][edge][1] + ctm.t[1];
609  b_1 = ctm.m[0] * ref_vert[mode][e->next_vert(edge)][0] + ctm.t[0];
610  b_2 = ctm.m[1] * ref_vert[mode][e->next_vert(edge)][1] + ctm.t[1];
611 
612  // values of nonpolynomial function in two vertices
613  double2 fa, fb;
614  calc_ref_map(e, nurbs, a_1, a_2, fa);
615  calc_ref_map(e, nurbs, b_1, b_2, fb);
616 
617  double2* pt = quad1d.get_points(mo1);
618  for (j = 0; j < np; j++) // over all integration points
619  {
620  double2 x, v;
621  double t = pt[j][0];
622  edge_coord(e, edge, t, x, v);
623  calc_ref_map(e, nurbs, x[0], x[1], fn[j]);
624 
625  for (k = 0; k < 2; k++)
626  fn[j][k] = fn[j][k] - (fa[k] + (t + 1)/2.0 * (fb[k] - fa[k]));
627  }
628 
629  double2* result = proj + e->get_nvert() + edge * (order - 1);
630  for (k = 0; k < 2; k++)
631  {
632  for (i = 0; i < ne; i++)
633  {
634  for (j = 0; j < np; j++)
635  {
636  double t = pt[j][0];
637  double fi = 0;
638  switch(i + 2)
639  {
640  case 0:
641  fi = l0(t);
642  break;
643  case 1:
644  fi = l1(t);
645  break;
646  case 2:
647  fi = l2(t);
648  break;
649  case 3:
650  fi = l3(t);
651  break;
652  case 4:
653  fi = l4(t);
654  break;
655  case 5:
656  fi = l5(t);
657  break;
658  case 6:
659  fi = l6(t);
660  break;
661  case 7:
662  fi = l7(t);
663  break;
664  case 8:
665  fi = l8(t);
666  break;
667  case 9:
668  fi = l9(t);
669  break;
670  case 10:
671  fi = l10(t);
672  break;
673  case 11:
674  fi = l11(t);
675  break;
676  }
677  rhside[k][i] += pt[j][1] * (fi * fn[j][k]);
678  }
679  }
680  // solve
681  cholsl(edge_proj_matrix, ne, edge_p, rhside[k], rhside[k]);
682  for (i = 0; i < ne; i++)
683  result[i][k] = rhside[k][i];
684  }
685  }
686 
688 
689  void CurvMap::old_projection(Element* e, int order, double2* proj, double* old[2], H1ShapesetJacobi* ref_map_shapeset, PrecalcShapeset* ref_map_pss)
690  {
691  int mo2 = quad2d.get_max_order(e->get_mode());
692  int np = quad2d.get_num_points(mo2, e->get_mode());
693 
694  for (unsigned int k = 0; k < e->get_nvert(); k++) // loop over vertices
695  {
696  // vertex basis functions in all integration points
697  double* vd;
698  int index_v = ref_map_shapeset->get_vertex_index(k, e->get_mode());
699  ref_map_pss->set_active_shape(index_v);
700  ref_map_pss->set_quad_order(mo2);
701  vd = ref_map_pss->get_fn_values();
702 
703  for (int m = 0; m < 2; m++) // part 0 or 1
704  for (int j = 0; j < np; j++)
705  old[m][j] += proj[k][m] * vd[j];
706 
707  for (int ii = 0; ii < order - 1; ii++)
708  {
709  // edge basis functions in all integration points
710  double* ed;
711  int index_e = ref_map_shapeset->get_edge_index(k, 0, ii + 2, e->get_mode());
712  ref_map_pss->set_active_shape(index_e);
713  ref_map_pss->set_quad_order(mo2);
714  ed = ref_map_pss->get_fn_values();
715 
716  for (int m = 0; m < 2; m++) //part 0 or 1
717  for (int j = 0; j < np; j++)
718  old[m][j] += proj[e->get_nvert() + k * (order-1) + ii][m] * ed[j];
719  }
720  }
721  }
722 
723  void CurvMap::calc_bubble_projection(Element* e, Nurbs** nurbs, int order, double2* proj, H1ShapesetJacobi* ref_map_shapeset, PrecalcShapeset* ref_map_pss)
724  {
725  ref_map_pss->set_active_element(e);
726 
727  int i, j, k;
728  int mo2 = quad2d.get_max_order(e->get_mode());
729  int np = quad2d.get_num_points(mo2, e->get_mode());
730  int qo = e->is_quad() ? H2D_MAKE_QUAD_ORDER(order, order) : order;
731  int nb = ref_map_shapeset->get_num_bubbles(qo, e->get_mode());
732 
733  double2* fn = new double2[np];
734  memset(fn, 0, np * sizeof(double2));
735 
736  double* rhside[2];
737  double* old[2];
738  for (i = 0; i < 2; i++)
739  {
740  rhside[i] = new double[nb];
741  old[i] = new double[np];
742  memset(rhside[i], 0, sizeof(double) * nb);
743  memset(old[i], 0, sizeof(double) * np);
744  }
745 
746  // compute known part of projection (vertex and edge part)
747  old_projection(e, order, proj, old, ref_map_shapeset, ref_map_pss);
748 
749  // fn values of both components of nonpolynomial function
750  double3* pt = quad2d.get_points(mo2, e->get_mode());
751  for (j = 0; j < np; j++) // over all integration points
752  {
753  double2 a;
754  a[0] = ctm.m[0] * pt[j][0] + ctm.t[0];
755  a[1] = ctm.m[1] * pt[j][1] + ctm.t[1];
756  calc_ref_map(e, nurbs, a[0], a[1], fn[j]);
757  }
758 
759  double2* result = proj + e->get_nvert() + e->get_nvert() * (order - 1);
760  for (k = 0; k < 2; k++)
761  {
762  for (i = 0; i < nb; i++) // loop over bubble basis functions
763  {
764  // bubble basis functions in all integration points
765  double *bfn;
766  int index_i = ref_map_shapeset->get_bubble_indices(qo, e->get_mode())[i];
767  ref_map_pss->set_active_shape(index_i);
768  ref_map_pss->set_quad_order(mo2);
769  bfn = ref_map_pss->get_fn_values();
770 
771  for (j = 0; j < np; j++) // over all integration points
772  rhside[k][i] += pt[j][2] * (bfn[j] * (fn[j][k] - old[k][j]));
773  }
774 
775  // solve
776  if(e->get_mode() == HERMES_MODE_TRIANGLE)
777  cholsl(bubble_proj_matrix_tri, nb, bubble_tri_p, rhside[k], rhside[k]);
778  else
779  cholsl(bubble_proj_matrix_quad, nb, bubble_quad_p, rhside[k], rhside[k]);
780 
781  for (i = 0; i < nb; i++)
782  result[i][k] = rhside[k][i];
783  }
784 
785  for (i = 0; i < 2; i++)
786  {
787  delete [] rhside[i];
788  delete [] old[i];
789  }
790  delete [] fn;
791  }
792 
794 
795  void CurvMap::ref_map_projection(Element* e, Nurbs** nurbs, int order, double2* proj, H1ShapesetJacobi* ref_map_shapeset, PrecalcShapeset* ref_map_pss)
796  {
797  // vertex part
798  for (unsigned int i = 0; i < e->get_nvert(); i++)
799  {
800  proj[i][0] = e->vn[i]->x;
801  proj[i][1] = e->vn[i]->y;
802  }
803 
804  if(e->cm->toplevel == false)
805  e = e->cm->parent;
806 
807  // edge part
808  for (int edge = 0; edge < (int)e->get_nvert(); edge++)
809  calc_edge_projection(e, edge, nurbs, order, proj, ref_map_shapeset, ref_map_pss);
810 
811  //bubble part
812  calc_bubble_projection(e, nurbs, order, proj, ref_map_shapeset, ref_map_pss);
813  }
814 
815  void CurvMap::update_refmap_coeffs(Element* e)
816  {
817  H1ShapesetJacobi ref_map_shapeset;
818  PrecalcShapeset ref_map_pss(&ref_map_shapeset);
819 
820  ref_map_pss.set_quad_2d(&quad2d);
821  ref_map_pss.set_active_element(e);
822 
823  // calculation of projection matrices
824  if(edge_proj_matrix == NULL)
825  precalculate_cholesky_projection_matrix_edge(&ref_map_shapeset, &ref_map_pss);
826  if(bubble_proj_matrix_tri == NULL && e->get_mode() == HERMES_MODE_TRIANGLE)
827  precalculate_cholesky_projection_matrices_bubble(&ref_map_shapeset, &ref_map_pss);
828  if(bubble_proj_matrix_quad == NULL && e->get_mode() == HERMES_MODE_QUAD)
829  precalculate_cholesky_projection_matrices_bubble(&ref_map_shapeset, &ref_map_pss);
830 
831  // allocate projection coefficients
832  int nv = e->get_nvert();
833  int ne = order - 1;
834  int qo = e->is_quad() ? H2D_MAKE_QUAD_ORDER(order, order) : order;
835  int nb = ref_map_shapeset.get_num_bubbles(qo, e->get_mode());
836  nc = nv + nv*ne + nb;
837  if(coeffs != NULL)
838  {
839  delete [] coeffs;
840  coeffs = NULL;
841  }
842  coeffs = new double2[nc];
843 
844  // WARNING: do not change the format of the array 'coeffs'. If it changes,
845  // RefMap::set_active_element() has to be changed too.
846 
847  Nurbs** nurbs;
848  if(toplevel == false)
849  {
850  ref_map_pss.set_active_element(e);
851  ref_map_pss.set_transform(part);
852  nurbs = parent->cm->nurbs;
853  }
854  else
855  {
856  ref_map_pss.reset_transform();
857  nurbs = e->cm->nurbs;
858  }
859  ctm = *(ref_map_pss.get_ctm());
860  ref_map_pss.reset_transform(); // fixme - do we need this?
861 
862  // calculation of new projection coefficients
863  ref_map_projection(e, nurbs, order, coeffs, &ref_map_shapeset, &ref_map_pss);
864  }
865 
866  void CurvMap::get_mid_edge_points(Element* e, double2* pt, int n)
867  {
868  Nurbs** nurbs = this->nurbs;
869  Transformable tran;
870  tran.set_active_element(e);
871 
872  if(toplevel == false)
873  {
874  tran.set_transform(part);
875  e = e->cm->parent;
876  nurbs = e->cm->nurbs;
877  }
878 
879  ctm = *(tran.get_ctm());
880  double xi_1, xi_2;
881  for (int i = 0; i < n; i++)
882  {
883  xi_1 = ctm.m[0] * pt[i][0] + ctm.t[0];
884  xi_2 = ctm.m[1] * pt[i][1] + ctm.t[1];
885  calc_ref_map(e, nurbs, xi_1, xi_2, pt[i]);
886  }
887  }
888 
889  void Nurbs::unref()
890  {
891  if(!--ref) // fixme: possible leak, we need ~Nurbs too
892  {
893  delete [] pt;
894  delete [] kv;
895  delete this;
896  }
897  }
898 
899  CurvMap::CurvMap(CurvMap* cm)
900  {
901  memcpy(this, cm, sizeof(CurvMap));
902  coeffs = new double2[nc];
903  memcpy(coeffs, cm->coeffs, sizeof(double2) * nc);
904 
905  if(toplevel)
906  for (int i = 0; i < 4; i++)
907  if(nurbs[i] != NULL)
908  nurbs[i]->ref++;
909  }
910 
911  CurvMap::~CurvMap()
912  {
913  if(coeffs != NULL)
914  {
915  delete [] coeffs;
916  coeffs = NULL;
917  }
918 
919  if(toplevel)
920  for (int i = 0; i < 4; i++)
921  if(nurbs[i] != NULL)
922  nurbs[i]->unref();
923  }
924  }
925 }