Hermes2D  2.0
shapeset.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 "global.h"
17 #include "shapeset.h"
18 #include "matrix.h"
19 
20 using namespace Hermes::Algebra::DenseMatrixOperations;
21 
22 namespace Hermes
23 {
24  namespace Hermes2D
25  {
26  double* Shapeset::calculate_constrained_edge_combination(int order, int part, int ori, ElementMode2D mode)
27  {
28  /*
29  "ebias" is the order of the first edge function, this has to be subtracted
30  from the order to get a reasonable numbering of the edge functions, starting
31  from 0. For H1 ebias is 2 and for Hcurl it is 0.
32  */
33  assert((order - ebias) >= 0);
34  assert(part >= 0);
35 
36  int i, j, n;
37 
38  // determine the interval of the edge
39  for (n = 2; n <= part; n <<= 1)
40  part -= n;
41 
42  double n2 = 2.0 / n;
43  double hi = -((double) part * n2 - 1.0);
44  double lo = -((double) (part + 1) * n2 - 1.0);
45 
46  int idx[16];
47  ori = ori ? 0 : 1;
48  for (i = 0; i <= order; i++)
49  idx[i] = get_edge_index(0, ori, i, mode);
50 
51  // function values at the endpoints (for subtracting of the linear part in H1)
52  double c = 1.0;
53  double f_lo = 0.0, f_hi = 0.0;
54  if(ebias == 2)
55  {
56  f_lo = get_value(0, idx[order], lo, -1.0, 0, mode);
57  f_hi = get_value(0, idx[order], hi, -1.0, 0, mode);
58  }
59  else
60  {
61  // this is for H(curl), where we need to multiply the constrained function
62  // so that its vectors are not shortened and fit the large constraining fn.
63  c = (hi - lo) / 2.0;
64  }
65 
66  // fill the matrix of the linear system
67  n = order + 1 - ebias;
68  int space_type = this->get_space_type();
69  int component = (space_type == HERMES_HDIV_SPACE)? 1 : 0;
70  double** a = new_matrix<double>(n, n);
71  double* b = new double[n];
72  for (i = 0; i < n; i++)
73  {
74  // chebyshev point
75  int o = (ebias == 0) ? order + 1 : order;
76  double p = cos((i + 1) * M_PI / o);
77  double r = (p + 1.0) * 0.5;
78  double s = 1.0 - r;
79 
80  // matrix row
81  for (j = 0; j < n; j++)
82  a[i][j] = get_value(0, idx[j + ebias], p, -1.0, component, mode);
83 
84  // rhs
85  b[i] = c * get_value(0, idx[order], lo*s + hi*r, -1.0, component, mode) - f_lo*s - f_hi*r;
86  }
87 
88  // solve the system
89  double d;
90  int* iperm = new int[n];
91  ludcmp(a, n, iperm, &d);
92  lubksb(a, n, iperm, b);
93 
94  // cleanup
95  delete [] iperm;
96  delete [] a;
97 
98  return b;
99  }
100 
101  double* Shapeset::get_constrained_edge_combination(int order, int part, int ori, int& nitems, ElementMode2D mode)
102  {
103  int index = 2*((max_order + 1 - ebias)*part + (order - ebias)) + ori;
104 
105  // allocate/reallocate the array if necessary
106  if(comb_table == NULL)
107  {
108  table_size = 1024;
109  while (table_size <= index) table_size *= 2;
110  comb_table = (double**) malloc(table_size * sizeof(double*));
111  memset(comb_table, 0, table_size * sizeof(double*));
112  }
113  else if(index >= table_size)
114  {
115  // adjust table_size to accommodate the required depth
116  int old_size = table_size;
117  while (index >= table_size) table_size *= 2;
118 
119  // reallocate the table
120  comb_table = (double**) realloc(comb_table, table_size * sizeof(double*));
121  memset(comb_table + old_size, 0, (table_size - old_size) * sizeof(double*));
122  }
123 
124  // do we have the required linear combination yet?
125  if(comb_table[index] == NULL)
126  {
127  // no, calculate it
128  comb_table[index] = calculate_constrained_edge_combination(order, part, ori, mode);
129  }
130 
131  nitems = order + 1 - ebias;
132  return comb_table[index];
133  }
134 
135  void Shapeset::free_constrained_edge_combinations()
136  {
137  if(comb_table != NULL)
138  {
139  for (int i = 0; i < table_size; i++)
140  if(comb_table[i] != NULL)
141  delete [] comb_table[i];
142 
143  free(comb_table);
144  comb_table = NULL;
145  }
146  }
147 
148  double Shapeset::get_constrained_value(int n, int index, double x, double y, int component, ElementMode2D mode)
149  {
150  index = -1 - index;
151 
152  int part = (unsigned) index >> 7;
153  int order = (index >> 3) & 15;
154  int edge = (index >> 1) & 3;
155  int ori = index & 1;
156 
157  int i, nc;
158  double sum, *comb = get_constrained_edge_combination(order, part, ori, nc, mode);
159 
160  sum = 0.0;
161  shape_fn_t* table = shape_table[n][mode][component];
162  for (i = 0; i < nc; i++)
163  sum += comb[i] * table[get_edge_index(edge, ori, i + ebias, mode)](x, y);
164 
165  return sum;
166  }
167 
168  Shapeset::~Shapeset() { free_constrained_edge_combinations(); }
169 
170  int Shapeset::get_max_order() const { return max_order; }
171 
172  int Shapeset::get_num_components() const { return num_components; }
173 
174  int Shapeset::get_vertex_index(int vertex, ElementMode2D mode) const
175  {
176  if(mode == HERMES_MODE_TRIANGLE)
177  assert(vertex < 3);
178  else
179  assert(vertex < 4);
180  return vertex_indices[mode][vertex];
181  }
182 
183  int Shapeset::get_edge_index(int edge, int ori, int order, ElementMode2D mode) const
184  {
185  if(mode == HERMES_MODE_TRIANGLE)
186  assert(edge < 3);
187  else
188  assert(edge < 4);
189  assert(order >= 0 && order <= max_order);
190  assert(ori == 0 || ori == 1);
191  return edge_indices[mode][edge][2*order + ori];
192  }
193 
194  int* Shapeset::get_bubble_indices(int order, ElementMode2D mode) const
195  {
196  assert(H2D_GET_H_ORDER(order) >= 0 && H2D_GET_H_ORDER(order) <= max_order);
197  assert(H2D_GET_V_ORDER(order) >= 0 && H2D_GET_V_ORDER(order) <= max_order);
198  int index = order;
199  if(mode == HERMES_MODE_QUAD) //tables of bubble indices are transposed
200  index = H2D_MAKE_QUAD_ORDER(H2D_GET_V_ORDER(order), H2D_GET_H_ORDER(order));
201  return bubble_indices[mode][index];
202  }
203 
204  int Shapeset::get_num_bubbles(int order, ElementMode2D mode) const
205  {
206  assert(H2D_GET_H_ORDER(order) >= 0 && H2D_GET_H_ORDER(order) <= max_order);
207  assert(H2D_GET_V_ORDER(order) >= 0 && H2D_GET_V_ORDER(order) <= max_order);
208  return bubble_count[mode][order];
209  }
210 
211  int Shapeset::get_constrained_edge_index(int edge, int order, int ori, int part, ElementMode2D mode) const
212  {
213  if(mode == HERMES_MODE_TRIANGLE)
214  assert(edge < 3);
215  else
216  assert(edge < 4);
217  assert(order >= 0 && order <= max_order);
218  assert(part >= 0);
219  assert(order <= H2D_ORDER_MASK);
220  return -1 - ((part << 7) + (order << 3) + (edge << 1) + ori);
221  }
222 
223  int Shapeset::get_order(int index, ElementMode2D mode) const
224  {
225  if(index >= 0) {
226  return index_to_order[mode][index];
227  }
228  else return ((-1 - index) >> 3) & 15;
229  }
230 
231  double Shapeset::get_value(int n, int index, double x, double y, int component, ElementMode2D mode)
232  {
233  if(index >= 0)
234  {
235  Shapeset::shape_fn_t** shape_expansion = shape_table[n][mode];
236  if(shape_expansion == NULL)
237  { // requested exansion (f, df/dx, df/dy, ddf/dxdx, ...) is not defined.
238  //just to keep the number of warnings low: warn just once about a given combinations of n, mode, and index.
239  static int warned_mode = -1, warned_index = -1, warned_n = 1;
240  this->warn_if(warned_mode != mode || warned_index != index || warned_n != n, "Requested undefined expansion %d (mode: %d) of a shape %d, returning 0", n, mode, index);
241  warned_mode = mode;
242  warned_index = index;
243  warned_n = n;
244  return 0.;
245  }
246  else
247  return shape_expansion[component][index](x, y);
248  }
249  else
250  return get_constrained_value(n, index, x, y, component, mode);
251  }
252 
253  double Shapeset::get_fn_value (int index, double x, double y, int component, ElementMode2D mode) { return get_value(0, index, x, y, component, mode); }
254  double Shapeset::get_dx_value (int index, double x, double y, int component, ElementMode2D mode) { return get_value(1, index, x, y, component, mode); }
255  double Shapeset::get_dy_value (int index, double x, double y, int component, ElementMode2D mode) { return get_value(2, index, x, y, component, mode); }
256  double Shapeset::get_dxx_value(int index, double x, double y, int component, ElementMode2D mode) { return get_value(3, index, x, y, component, mode); }
257  double Shapeset::get_dyy_value(int index, double x, double y, int component, ElementMode2D mode) { return get_value(4, index, x, y, component, mode); }
258  double Shapeset::get_dxy_value(int index, double x, double y, int component, ElementMode2D mode) { return get_value(5, index, x, y, component, mode); }
259 
261  double2* Shapeset::get_ref_vertex(int vertex, ElementMode2D mode)
262  {
263  return &ref_vert[mode][vertex];
264  }
265  }
266 }