Hermes2D  2.0
h1_proj_based_selector.cpp
1 #include "global.h"
2 #include "matrix.h"
3 #include "solution.h"
4 #include "shapeset/shapeset_h1_all.h"
5 #include "element_to_refine.h"
6 #include "h1_proj_based_selector.h"
7 namespace Hermes
8 {
9  namespace Hermes2D
10  {
11  namespace RefinementSelectors
12  {
13  template<typename Scalar>
14  const int H1ProjBasedSelector<Scalar>::H2DRS_MAX_H1_ORDER = H2DRS_MAX_ORDER;
15 
16  template<typename Scalar>
17  H1ProjBasedSelector<Scalar>::H1ProjBasedSelector(CandList cand_list, double conv_exp, int max_order, H1Shapeset* user_shapeset)
18  : ProjBasedSelector<Scalar>(cand_list, conv_exp, max_order, user_shapeset == NULL ? new H1Shapeset() : user_shapeset, typename OptimumSelector<Scalar>::Range(1, 1), typename OptimumSelector<Scalar>::Range(2, H2DRS_MAX_H1_ORDER)), user_shapeset(user_shapeset == NULL ? false : true)
19  {
20  if(user_shapeset != NULL)
21  {
22  this->warn("Warn: The user shapeset provided for the selector has to have a correct copy constructor implemented.");
23  this->warn("Warning: The functionality for cloning user shapeset is to be implemented yet.");
24  }
25  }
26 
27  template<typename Scalar>
28  H1ProjBasedSelector<Scalar>::~H1ProjBasedSelector()
29  {
30  if(!this->user_shapeset)
31  delete this->shapeset;
32  }
33 
34  template<typename Scalar>
35  Selector<Scalar>* H1ProjBasedSelector<Scalar>::clone()
36  {
37  H1ProjBasedSelector<Scalar>* newSelector = new H1ProjBasedSelector(this->cand_list, this->conv_exp, this->max_order, (H1Shapeset*)this->shapeset);
38  newSelector->set_error_weights(this->error_weight_h, this->error_weight_p, this->error_weight_aniso);
39  newSelector->isAClone = true;
40  return newSelector;
41  }
42 
43  template<typename Scalar>
44  void H1ProjBasedSelector<Scalar>::set_current_order_range(Element* element)
45  {
46  this->current_max_order = this->max_order;
47  int max_element_order = (20 - element->iro_cache)/2 - 1;
48  if(this->current_max_order == H2DRS_DEFAULT_ORDER)
49  this->current_max_order = max_element_order; // default
50  else
51  this->current_max_order = std::min(this->current_max_order, max_element_order); // user specified
52  this->current_min_order = 1;
53  }
54 
55  template<typename Scalar>
56  void H1ProjBasedSelector<Scalar>::precalc_shapes(const double3* gip_points, const int num_gip_points, const Trf* trfs, const int num_noni_trfs, const Hermes::vector<typename OptimumSelector<Scalar>::ShapeInx>& shapes, const int max_shape_inx, typename ProjBasedSelector<Scalar>::TrfShape& svals, ElementMode2D mode)
57  {
58  //for all transformations
59  bool done = false;
60  int inx_trf = 0;
61  while (!done && inx_trf < H2D_TRF_NUM)
62  {
63  //prepare data for processing
64  const Trf& trf = trfs[inx_trf];
65  Hermes::vector<typename ProjBasedSelector<Scalar>::TrfShapeExp>& trf_svals = svals[inx_trf];
66 
67  //allocate
68  trf_svals.resize(max_shape_inx + 1);
69 
70  //for all shapes
71  const int num_shapes = (int)shapes.size();
72  for(int i = 0; i < num_shapes; i++)
73  {
74  int inx_shape = shapes[i].inx;
75  typename ProjBasedSelector<Scalar>::TrfShapeExp& shape_exp = trf_svals[inx_shape];
76 
77  //allocate
78  shape_exp.allocate(H2D_H1FE_NUM, num_gip_points);
79 
80  //for all GIP points
81  for(int k = 0; k < num_gip_points; k++)
82  {
83  //transform coordinates
84  double ref_x = gip_points[k][H2D_GIP2D_X] * trf.m[0] + trf.t[0];
85  double ref_y = gip_points[k][H2D_GIP2D_Y] * trf.m[1] + trf.t[1];
86 
87  //for all expansions: retrieve values
88  shape_exp[H2D_H1FE_VALUE][k] = this->shapeset->get_fn_value(inx_shape, ref_x, ref_y, 0, mode);
89  shape_exp[H2D_H1FE_DX][k] = this->shapeset->get_dx_value(inx_shape, ref_x, ref_y, 0, mode);
90  shape_exp[H2D_H1FE_DY][k] = this->shapeset->get_dy_value(inx_shape, ref_x, ref_y, 0, mode);
91  }
92  }
93 
94  //move to the next transformation
95  if(inx_trf == H2D_TRF_IDENTITY)
96  done = true;
97  else
98  {
99  inx_trf++;
100  if(inx_trf >= num_noni_trfs) //if all transformations were processed, move to the identity transformation
101  inx_trf = H2D_TRF_IDENTITY;
102  }
103  }
104  if(!done)
105  throw Exceptions::Exception("All transformation processed but identity transformation not found."); //identity transformation has to be the last transformation
106  }
107 
108  template<typename Scalar>
109  void H1ProjBasedSelector<Scalar>::precalc_ortho_shapes(const double3* gip_points, const int num_gip_points, const Trf* trfs, const int num_noni_trfs, const Hermes::vector<typename OptimumSelector<Scalar>::ShapeInx>& shapes, const int max_shape_inx, typename ProjBasedSelector<Scalar>::TrfShape& svals, ElementMode2D mode)
110  {
111  //calculate values
112  precalc_shapes(gip_points, num_gip_points, trfs, num_noni_trfs, shapes, max_shape_inx, svals, mode);
113 
114  //calculate orthonormal basis
115  const int num_shapes = (int)shapes.size();
116  for(int i = 0; i < num_shapes; i++)
117  {
118  const int inx_shape_i = shapes[i].inx;
119 
120  //orthogonalize
121  for(int j = 0; j < i; j++)
122  {
123  const int inx_shape_j = shapes[j].inx;
124 
125  //calculate product of non-transformed functions
126  double product = 0.0;
127  for(int k = 0; k < num_gip_points; k++)
128  {
129  double sum = 0.0;
130  sum += svals[H2D_TRF_IDENTITY][inx_shape_i][H2D_H1FE_VALUE][k] * svals[H2D_TRF_IDENTITY][inx_shape_j][H2D_H1FE_VALUE][k];
131  sum += svals[H2D_TRF_IDENTITY][inx_shape_i][H2D_H1FE_DX][k] * svals[H2D_TRF_IDENTITY][inx_shape_j][H2D_H1FE_DX][k];
132  sum += svals[H2D_TRF_IDENTITY][inx_shape_i][H2D_H1FE_DY][k] * svals[H2D_TRF_IDENTITY][inx_shape_j][H2D_H1FE_DY][k];
133  product += gip_points[k][H2D_GIP2D_W] * sum;
134  }
135 
136  //for all transformations
137  int inx_trf = 0;
138  bool done = false;
139  while (!done && inx_trf < H2D_TRF_NUM)
140  {
141  //for all integration points
142  for(int k = 0; k < num_gip_points; k++)
143  {
144  svals[inx_trf][inx_shape_i][H2D_H1FE_VALUE][k] -= product * svals[inx_trf][inx_shape_j][H2D_H1FE_VALUE][k];
145  svals[inx_trf][inx_shape_i][H2D_H1FE_DX][k] -= product * svals[inx_trf][inx_shape_j][H2D_H1FE_DX][k];
146  svals[inx_trf][inx_shape_i][H2D_H1FE_DY][k] -= product * svals[inx_trf][inx_shape_j][H2D_H1FE_DY][k];
147  }
148 
149  //move to the next transformation
150  if(inx_trf == H2D_TRF_IDENTITY)
151  done = true;
152  else
153  {
154  inx_trf++;
155  if(inx_trf >= num_noni_trfs) //if all transformations were processed, move to the identity transformation
156  inx_trf = H2D_TRF_IDENTITY;
157  }
158  }
159  if(!done)
160  throw Exceptions::Exception("All transformation processed but identity transformation not found."); //identity transformation has to be the last transformation
161  }
162 
163  //normalize
164  //calculate norm
165  double norm_squared = 0.0;
166  for(int k = 0; k < num_gip_points; k++)
167  {
168  double sum = 0.0;
169  sum += sqr(svals[H2D_TRF_IDENTITY][inx_shape_i][H2D_H1FE_VALUE][k]);
170  sum += sqr(svals[H2D_TRF_IDENTITY][inx_shape_i][H2D_H1FE_DX][k]);
171  sum += sqr(svals[H2D_TRF_IDENTITY][inx_shape_i][H2D_H1FE_DY][k]);
172  norm_squared += gip_points[k][H2D_GIP2D_W] * sum;
173  }
174  double norm = sqrt(norm_squared);
175  if(!finite(1/norm))
176  throw Exceptions::Exception("Norm (%g) is almost zero.", norm);
177 
178  //for all transformations: normalize
179  int inx_trf = 0;
180  bool done = false;
181  while (!done && inx_trf < H2D_TRF_NUM)
182  {
183  //for all integration points
184  for(int k = 0; k < num_gip_points; k++)
185  {
186  svals[inx_trf][inx_shape_i][H2D_H1FE_VALUE][k] /= norm;
187  svals[inx_trf][inx_shape_i][H2D_H1FE_DX][k] /= norm;
188  svals[inx_trf][inx_shape_i][H2D_H1FE_DY][k] /= norm;
189  }
190 
191  //move to the next transformation
192  if(inx_trf == H2D_TRF_IDENTITY)
193  done = true;
194  else
195  {
196  inx_trf++;
197  if(inx_trf >= num_noni_trfs) //if all transformations were processed, move to the identity transformation
198  inx_trf = H2D_TRF_IDENTITY;
199  }
200  }
201  if(!done)
202  throw Exceptions::Exception("All transformation processed but identity transformation not found."); //identity transformation has to be the last transformation
203  }
204  }
205 
206  template<typename Scalar>
207  Scalar** H1ProjBasedSelector<Scalar>::precalc_ref_solution(int inx_son, Solution<Scalar>* rsln, Element* element, int intr_gip_order)
208  {
209  //fill with values
210  Scalar** rvals_son = precalc_rvals[inx_son];
211  rvals_son[H2D_H1FE_VALUE] = rsln->get_fn_values(0);
212  rvals_son[H2D_H1FE_DX] = rsln->get_dx_values(0);
213  rvals_son[H2D_H1FE_DY] = rsln->get_dy_values(0);
214 
215  return rvals_son;
216  }
217 
218  template<typename Scalar>
219  double** H1ProjBasedSelector<Scalar>::build_projection_matrix(double3* gip_points, int num_gip_points,
220  const int* shape_inx, const int num_shapes, ElementMode2D mode)
221  {
222  //allocate
223  double** matrix = new_matrix<double>(num_shapes, num_shapes);
224 
225  //calculate products
226  int inx_row = 0;
227  for(int i = 0; i < num_shapes; i++, inx_row += num_shapes)
228  {
229  double* matrix_row = matrix[i];
230  int shape0_inx = shape_inx[i];
231  for(int k = 0; k < num_shapes; k++)
232  {
233  int shape1_inx = shape_inx[k];
234 
235  double value = 0.0;
236  for(int j = 0; j < num_gip_points; j++)
237  {
238  double gip_x = gip_points[j][H2D_GIP2D_X], gip_y = gip_points[j][H2D_GIP2D_Y];
239  double value0 = this->shapeset->get_value(H2D_FEI_VALUE, shape0_inx, gip_x, gip_y, 0, mode);
240  double value1 = this->shapeset->get_value(H2D_FEI_VALUE, shape1_inx, gip_x, gip_y, 0, mode);
241  double dx0 = this->shapeset->get_value(H2D_FEI_DX, shape0_inx, gip_x, gip_y, 0, mode);
242  double dx1 = this->shapeset->get_value(H2D_FEI_DX, shape1_inx, gip_x, gip_y, 0, mode);
243  double dy0 = this->shapeset->get_value(H2D_FEI_DY, shape0_inx, gip_x, gip_y, 0, mode);
244  double dy1 = this->shapeset->get_value(H2D_FEI_DY, shape1_inx, gip_x, gip_y, 0, mode);
245 
246  value += gip_points[j][H2D_GIP2D_W] * (value0*value1 + dx0*dx1 + dy0*dy1);
247  }
248 
249  matrix_row[k] = value;
250  }
251  }
252 
253  return matrix;
254  }
255 
256  template<typename Scalar>
257  Scalar H1ProjBasedSelector<Scalar>::evaluate_rhs_subdomain(Element* sub_elem, const typename ProjBasedSelector<Scalar>::ElemGIP& sub_gip, const typename ProjBasedSelector<Scalar>::ElemSubTrf& sub_trf, const typename ProjBasedSelector<Scalar>::ElemSubShapeFunc& sub_shape)
258  {
259  Scalar total_value = 0;
260  for(int gip_inx = 0; gip_inx < sub_gip.num_gip_points; gip_inx++)
261  {
262  double3 &gip_pt = sub_gip.gip_points[gip_inx];
263 
264  //get value of a shape function
265  double shape_value[H2D_H1FE_NUM] = {0, 0, 0};
266  shape_value[H2D_H1FE_VALUE] = sub_shape.svals[H2D_H1FE_VALUE][gip_inx];
267  shape_value[H2D_H1FE_DX] = sub_shape.svals[H2D_H1FE_DX][gip_inx];
268  shape_value[H2D_H1FE_DY] = sub_shape.svals[H2D_H1FE_DY][gip_inx];
269 
271  //double ref_x = gip_pt[H2D_GIP2D_X] * sub_trf.trf->m[0] + sub_trf.trf->t[0];
272  //double ref_y = gip_pt[H2D_GIP2D_Y] * sub_trf.trf->m[1] + sub_trf.trf->t[1];
273  //Scalar shape_value[H2D_H1FE_NUM] = {0, 0, 0};
274  //shape_value[H2D_H1FE_VALUE] = shapeset->get_fn_value(sub_shape.inx, ref_x, ref_y, 0);
275  //shape_value[H2D_H1FE_DX] = shapeset->get_dx_value(sub_shape.inx, ref_x, ref_y, 0);
276  //shape_value[H2D_H1FE_DY] = shapeset->get_dy_value(sub_shape.inx, ref_x, ref_y, 0);
277  //error_if(std::abs(shape_value[H2D_H1FE_VALUE] - shape_valueA[H2D_H1FE_VALUE]) > 1E-15
278  // || std::abs(shape_value[H2D_H1FE_DX] - shape_valueA[H2D_H1FE_DX]) > 1E-15
279  // || std::abs(shape_value[H2D_H1FE_DY] - shape_valueA[H2D_H1FE_DY]) > 1E-15, "A1");
281 
282  //get value of ref. solution
283  Scalar ref_value[H2D_H1FE_NUM];
284  ref_value[H2D_H1FE_VALUE] = sub_gip.rvals[H2D_H1FE_VALUE][gip_inx];
285  ref_value[H2D_H1FE_DX] = sub_trf.coef_mx * sub_gip.rvals[H2D_H1FE_DX][gip_inx];
286  ref_value[H2D_H1FE_DY] = sub_trf.coef_my * sub_gip.rvals[H2D_H1FE_DY][gip_inx];
287 
288  //evaluate a right-hand value
289  Scalar value = (shape_value[H2D_H1FE_VALUE] * ref_value[H2D_H1FE_VALUE])
290  + (shape_value[H2D_H1FE_DX] * ref_value[H2D_H1FE_DX])
291  + (shape_value[H2D_H1FE_DY] * ref_value[H2D_H1FE_DY]);
292 
293  total_value += gip_pt[H2D_GIP2D_W] * value;
294  }
295  return total_value;
296  }
297 
298  template<typename Scalar>
299  double H1ProjBasedSelector<Scalar>::evaluate_error_squared_subdomain(Element* sub_elem, const typename ProjBasedSelector<Scalar>::ElemGIP& sub_gip, const typename ProjBasedSelector<Scalar>::ElemSubTrf& sub_trf, const typename ProjBasedSelector<Scalar>::ElemProj& elem_proj)
300  {
301  double total_error_squared = 0;
302  for(int gip_inx = 0; gip_inx < sub_gip.num_gip_points; gip_inx++)
303  {
304  double3 &gip_pt = sub_gip.gip_points[gip_inx];
305 
306  //calculate value of projected solution
307  Scalar proj_value[H2D_H1FE_NUM] = {0, 0, 0};
308  for(int i = 0; i < elem_proj.num_shapes; i++)
309  {
310  int shape_inx = elem_proj.shape_inxs[i];
311  proj_value[H2D_H1FE_VALUE] += elem_proj.shape_coeffs[i] * elem_proj.svals[shape_inx][H2D_H1FE_VALUE][gip_inx];
312  proj_value[H2D_H1FE_DX] += elem_proj.shape_coeffs[i] * elem_proj.svals[shape_inx][H2D_H1FE_DX][gip_inx];
313  proj_value[H2D_H1FE_DY] += elem_proj.shape_coeffs[i] * elem_proj.svals[shape_inx][H2D_H1FE_DY][gip_inx];
314  }
315 
317  //double ref_x = gip_pt[H2D_GIP2D_X] * sub_trf.trf->m[0] + sub_trf.trf->t[0];
318  //double ref_y = gip_pt[H2D_GIP2D_Y] * sub_trf.trf->m[1] + sub_trf.trf->t[1];
319  //Scalar proj_valueA[H2D_H1FE_NUM] = {0, 0, 0};
320  //for(int i = 0; i < elem_proj.num_shapes; i++)
321  {
322  // int shape_inx = elem_proj.shape_inxs[i];
323  // proj_valueA[H2D_H1FE_VALUE] += elem_proj.shape_coeffs[i] * shapeset->get_fn_value(shape_inx, ref_x, ref_y, 0);
324  // proj_valueA[H2D_H1FE_DX] += elem_proj.shape_coeffs[i] * shapeset->get_dx_value(shape_inx, ref_x, ref_y, 0);
325  // proj_valueA[H2D_H1FE_DY] += elem_proj.shape_coeffs[i] * shapeset->get_dy_value(shape_inx, ref_x, ref_y, 0);
326  //}
327  //error_if(std::abs(proj_value[H2D_H1FE_VALUE] - proj_valueA[H2D_H1FE_VALUE]) > 1E-15
328  // || std::abs(proj_value[H2D_H1FE_DX] - proj_valueA[H2D_H1FE_DX]) > 1E-15
329  // || std::abs(proj_value[H2D_H1FE_DY] - proj_valueA[H2D_H1FE_DY]) > 1E-15, "A2");
331 
332  //get value of ref. solution
333  Scalar ref_value[3];
334  ref_value[H2D_H1FE_VALUE] = sub_gip.rvals[H2D_H1FE_VALUE][gip_inx];
335  ref_value[H2D_H1FE_DX] = sub_trf.coef_mx * sub_gip.rvals[H2D_H1FE_DX][gip_inx];
336  ref_value[H2D_H1FE_DY] = sub_trf.coef_my * sub_gip.rvals[H2D_H1FE_DY][gip_inx];
337 
338  //evaluate error
339  double error_squared = sqr(proj_value[H2D_H1FE_VALUE] - ref_value[H2D_H1FE_VALUE])
340  + sqr(proj_value[H2D_H1FE_DX] - ref_value[H2D_H1FE_DX])
341  + sqr(proj_value[H2D_H1FE_DY] - ref_value[H2D_H1FE_DY]);
342 
343  total_error_squared += gip_pt[H2D_GIP2D_W] * error_squared;
344  }
345  }
346  return total_error_squared;
347  }
348  template class HERMES_API H1ProjBasedSelector<double>;
349  template class HERMES_API H1ProjBasedSelector<std::complex<double> >;
350  }
351  }
352 }