Hermes2D  3.0
hcurl_proj_based_selector.cpp
1 #include "matrix.h"
2 #include "hcurl_proj_based_selector.h"
3 
4 namespace Hermes
5 {
6  namespace Hermes2D
7  {
8  namespace RefinementSelectors
9  {
10  template<typename Scalar>
11  const int HcurlProjBasedSelector<Scalar>::H2DRS_MAX_HCURL_ORDER = 6;
12 
13  template<typename Scalar>
15  : ProjBasedSelector<Scalar>(cand_list, max_order, user_shapeset == nullptr ? new HcurlShapeset() : user_shapeset, Range(), Range(0, H2DRS_MAX_HCURL_ORDER))
16  {
17  if (user_shapeset != nullptr)
18  {
19  this->warn("Warning: The user shapeset provided for the selector has to have a correct copy constructor implemented.");
20  this->warn("Warning: The functionality for cloning user shapeset is to be implemented yet.");
21  }
22  }
23 
24  template<typename Scalar>
26  {
27  }
28 
29  template<typename Scalar>
30  void HcurlProjBasedSelector<Scalar>::get_current_order_range(Element* element, int& min_order_, int& max_order_)
31  {
32  int max_element_order = (20 - element->iro_cache) / 2 - 1;
33  if (this->max_order == H2DRS_DEFAULT_ORDER)
34  max_order_ = max_element_order;
35  else
36  max_order_ = std::min(this->max_order, max_element_order);
37  min_order_ = 0;
38  }
39 
40  template<typename Scalar>
41  void HcurlProjBasedSelector<Scalar>::precalc_shapes(const double3* gip_points, const int num_gip_points, const Trf* trfs, const int num_noni_trfs, const std::vector<typename OptimumSelector<Scalar>::ShapeInx>& shapes, const int max_shape_inx, typename ProjBasedSelector<Scalar>::TrfShape& svals, ElementMode2D mode)
42  {
43  //for all transformations
44  bool done = false;
45  int inx_trf = 0;
46  while (!done && inx_trf < H2D_TRF_NUM)
47  {
48  //prepare data for processing
49  const Trf& trf = trfs[inx_trf];
50  std::vector<typename ProjBasedSelector<Scalar>::TrfShapeExp>& trf_svals = svals[inx_trf];
51 
52  //allocate
53  trf_svals.resize(max_shape_inx + 1);
54 
55  //for all shapes
56  const int num_shapes = (int)shapes.size();
57  for (int i = 0; i < num_shapes; i++)
58  {
59  int inx_shape = shapes[i].inx;
60  typename ProjBasedSelector<Scalar>::TrfShapeExp& shape_exp = trf_svals[inx_shape];
61 
62  //allocate
63  shape_exp.allocate(H2D_HCFE_NUM, num_gip_points);
64 
65  //for all GIP points
66  for (int k = 0; k < num_gip_points; k++)
67  {
68  //transform coordinates
69  double ref_x = gip_points[k][H2D_GIP2D_X] * trf.m[0] + trf.t[0];
70  double ref_y = gip_points[k][H2D_GIP2D_Y] * trf.m[1] + trf.t[1];
71 
72  //for all expansions: retrieve values
73  shape_exp[H2D_HCFE_VALUE0][k] = this->shapeset->get_fn_value(inx_shape, ref_x, ref_y, 0, mode);
74  shape_exp[H2D_HCFE_VALUE1][k] = this->shapeset->get_fn_value(inx_shape, ref_x, ref_y, 1, mode);
75  shape_exp[H2D_HCFE_CURL][k] = this->shapeset->get_dx_value(inx_shape, ref_x, ref_y, 1, mode) - this->shapeset->get_dy_value(inx_shape, ref_x, ref_y, 0, mode);
76  }
77  }
78 
79  //move to the next transformation
80  if (inx_trf == H2D_TRF_IDENTITY)
81  done = true;
82  else
83  {
84  inx_trf++;
85  if (inx_trf >= num_noni_trfs) //if all transformations were processed, move to the identity transformation
86  inx_trf = H2D_TRF_IDENTITY;
87  }
88  }
89  if (!done)
90  //identity transformation has to be the last transformation
91  throw Exceptions::Exception("All transformation processed but identity transformation not found.");
92  }
93 
94  template<typename Scalar>
95  void HcurlProjBasedSelector<Scalar>::precalc_ortho_shapes(const double3* gip_points, const int num_gip_points, const Trf* trfs, const int num_noni_trfs, const std::vector<typename OptimumSelector<Scalar>::ShapeInx>& shapes, const int max_shape_inx, typename ProjBasedSelector<Scalar>::TrfShape& svals, ElementMode2D mode)
96  {
97  //calculate values
98  precalc_shapes(gip_points, num_gip_points, trfs, num_noni_trfs, shapes, max_shape_inx, svals, mode);
99 
100  //calculate orthonormal basis
101  const int num_shapes = (int)shapes.size();
102  for (int i = 0; i < num_shapes; i++)
103  {
104  const int inx_shape_i = shapes[i].inx;
105 
106  //orthogonalize
107  for (int j = 0; j < i; j++)
108  {
109  const int inx_shape_j = shapes[j].inx;
110 
111  //calculate product of non-transformed functions
112  double product = 0.0;
113  for (int k = 0; k < num_gip_points; k++)
114  {
115  double sum = 0.0;
116  sum += svals[H2D_TRF_IDENTITY][inx_shape_i][H2D_HCFE_VALUE0][k] * svals[H2D_TRF_IDENTITY][inx_shape_j][H2D_HCFE_VALUE0][k];
117  sum += svals[H2D_TRF_IDENTITY][inx_shape_i][H2D_HCFE_VALUE1][k] * svals[H2D_TRF_IDENTITY][inx_shape_j][H2D_HCFE_VALUE1][k];
118  sum += svals[H2D_TRF_IDENTITY][inx_shape_i][H2D_HCFE_CURL][k] * svals[H2D_TRF_IDENTITY][inx_shape_j][H2D_HCFE_CURL][k];
119  product += gip_points[k][H2D_GIP2D_W] * sum;
120  }
121 
122  //for all transformations
123  int inx_trf = 0;
124  bool done = false;
125  while (!done && inx_trf < H2D_TRF_NUM)
126  {
127  //for all integration points
128  for (int k = 0; k < num_gip_points; k++)
129  {
130  svals[inx_trf][inx_shape_i][H2D_HCFE_VALUE0][k] -= product * svals[inx_trf][inx_shape_j][H2D_HCFE_VALUE0][k];
131  svals[inx_trf][inx_shape_i][H2D_HCFE_VALUE1][k] -= product * svals[inx_trf][inx_shape_j][H2D_HCFE_VALUE1][k];
132  svals[inx_trf][inx_shape_i][H2D_HCFE_CURL][k] -= product * svals[inx_trf][inx_shape_j][H2D_HCFE_CURL][k];
133  }
134 
135  //move to the next transformation
136  if (inx_trf == H2D_TRF_IDENTITY)
137  done = true;
138  else
139  {
140  inx_trf++;
141  if (inx_trf >= num_noni_trfs) //if all transformations were processed, move to the identity transformation
142  inx_trf = H2D_TRF_IDENTITY;
143  }
144  }
145  if (!done)
146  //identity transformation has to be the last transformation
147  throw Exceptions::Exception("All transformation processed but identity transformation not found.");
148  }
149 
150  //normalize
151  //calculate norm
152  double norm_squared = 0.0;
153  for (int k = 0; k < num_gip_points; k++)
154  {
155  double sum = 0.0;
156  sum += sqr(svals[H2D_TRF_IDENTITY][inx_shape_i][H2D_HCFE_VALUE0][k]);
157  sum += sqr(svals[H2D_TRF_IDENTITY][inx_shape_i][H2D_HCFE_VALUE1][k]);
158  sum += sqr(svals[H2D_TRF_IDENTITY][inx_shape_i][H2D_HCFE_CURL][k]);
159  norm_squared += gip_points[k][H2D_GIP2D_W] * sum;
160  }
161  double norm = sqrt(norm_squared);
162  if (!finite(1 / norm))
163  throw Exceptions::Exception("Norm (%g) is almost zero.", norm);
164 
165  //for all transformations: normalize
166  int inx_trf = 0;
167  bool done = false;
168  while (!done && inx_trf < H2D_TRF_NUM)
169  {
170  //for all integration points
171  for (int k = 0; k < num_gip_points; k++)
172  {
173  svals[inx_trf][inx_shape_i][H2D_HCFE_VALUE0][k] /= norm;
174  svals[inx_trf][inx_shape_i][H2D_HCFE_VALUE1][k] /= norm;
175  svals[inx_trf][inx_shape_i][H2D_HCFE_CURL][k] /= norm;
176  }
177 
178  //move to the next transformation
179  if (inx_trf == H2D_TRF_IDENTITY)
180  done = true;
181  else
182  {
183  inx_trf++;
184  if (inx_trf >= num_noni_trfs) //if all transformations were processed, move to the identity transformation
185  inx_trf = H2D_TRF_IDENTITY;
186  }
187  }
188  if (!done)
189  //identity transformation has to be the last transformation
190  throw Exceptions::Exception("All transformation processed but identity transformation not found.");
191  }
192  }
193 
194  template<typename Scalar>
195  void HcurlProjBasedSelector<Scalar>::precalc_ref_solution(int inx_son, MeshFunction<Scalar>* rsln, Element* element, int intr_gip_order, Scalar* rval[H2D_MAX_ELEMENT_SONS][MAX_NUMBER_FUNCTION_VALUES_FOR_SELECTORS])
196  {
197  const int num_gip = rsln->get_quad_2d()->get_num_points(intr_gip_order, rsln->get_active_element()->get_mode());
198 
199  rval[inx_son][H2D_HCFE_VALUE0] = malloc_with_check<Scalar>(num_gip);
200  rval[inx_son][H2D_HCFE_VALUE1] = malloc_with_check<Scalar>(num_gip);
201  rval[inx_son][H2D_HCFE_CURL] = malloc_with_check<Scalar>(num_gip);
202 
203  memcpy(rval[inx_son][H2D_HCFE_VALUE0], rsln->get_fn_values(0), num_gip * sizeof(Scalar));
204  memcpy(rval[inx_son][H2D_HCFE_VALUE1], rsln->get_fn_values(1), num_gip * sizeof(Scalar));
205 
206  // prepare for curl
207  const Scalar* d1dx = rsln->get_dx_values(1);
208  const Scalar* d0dy = rsln->get_dy_values(0);
209  for (int i = 0; i < num_gip; i++)
210  rval[inx_son][H2D_HCFE_CURL][i] = d1dx[i] - d0dy[i];
211  }
212 
213  template<typename Scalar>
214  void HcurlProjBasedSelector<Scalar>::free_ref_solution_data(int inx_son, Scalar* rval[H2D_MAX_ELEMENT_SONS][MAX_NUMBER_FUNCTION_VALUES_FOR_SELECTORS])
215  {
216  free_with_check(rval[inx_son][H2D_HCFE_VALUE0]);
217  free_with_check(rval[inx_son][H2D_HCFE_VALUE1]);
218  free_with_check(rval[inx_son][H2D_HCFE_CURL]);
219  }
220 
221  template<typename Scalar>
222  double** HcurlProjBasedSelector<Scalar>::build_projection_matrix(double3* gip_points, int num_gip_points,
223  const int* shape_inx, const int num_shapes, ElementMode2D mode)
224  {
225  //allocate
226  double** matrix = new_matrix<double>(num_shapes, num_shapes);
227 
228  //calculate products
229  int inx_row = 0;
230  for (int i = 0; i < num_shapes; i++, inx_row += num_shapes)
231  {
232  double* matrix_row = matrix[i];
233  int shape0_inx = shape_inx[i];
234  for (int k = 0; k < num_shapes; k++)
235  {
236  int shape1_inx = shape_inx[k];
237 
238  double value = 0.0;
239  for (int j = 0; j < num_gip_points; j++)
240  {
241  double gip_x = gip_points[j][H2D_GIP2D_X], gip_y = gip_points[j][H2D_GIP2D_Y];
242  double value0[2] = { this->shapeset->get_value(H2D_FEI_VALUE, shape0_inx, gip_x, gip_y, 0, mode), this->shapeset->get_value(H2D_FEI_VALUE, shape0_inx, gip_x, gip_y, 1, mode) };
243  double value1[2] = { this->shapeset->get_value(H2D_FEI_VALUE, shape1_inx, gip_x, gip_y, 0, mode), this->shapeset->get_value(H2D_FEI_VALUE, shape1_inx, gip_x, gip_y, 1, mode) };
244  double d1dx0 = this->shapeset->get_value(H2D_FEI_DX, shape0_inx, gip_x, gip_y, 1, mode);
245  double d1dx1 = this->shapeset->get_value(H2D_FEI_DX, shape1_inx, gip_x, gip_y, 1, mode);
246  double d0dy0 = this->shapeset->get_value(H2D_FEI_DY, shape0_inx, gip_x, gip_y, 0, mode);
247  double d0dy1 = this->shapeset->get_value(H2D_FEI_DY, shape1_inx, gip_x, gip_y, 0, mode);
248  double curl0 = d1dx0 - d0dy0;
249  double curl1 = d1dx1 - d0dy1;
250 
251  value += gip_points[j][H2D_GIP2D_W] * (value0[0] * value1[0] + value0[1] * value1[1] + curl0*curl1);
252  }
253 
254  matrix_row[k] = value;
255  }
256  }
257 
258  return matrix;
259  }
260 
261  template<typename Scalar>
262  Scalar HcurlProjBasedSelector<Scalar>::evaluate_rhs_subdomain(Element* sub_elem, const typename ProjBasedSelector<Scalar>::ElemGIP& sub_gip, int son, const typename ProjBasedSelector<Scalar>::ElemSubTrf& sub_trf, const typename ProjBasedSelector<Scalar>::ElemSubShapeFunc& sub_shape, Scalar* rval[H2D_MAX_ELEMENT_SONS][MAX_NUMBER_FUNCTION_VALUES_FOR_SELECTORS])
263  {
264  double coef_curl = std::abs(sub_trf.coef_mx * sub_trf.coef_my);
265  Scalar total_value = 0;
266  for (int gip_inx = 0; gip_inx < sub_gip.num_gip_points; gip_inx++)
267  {
268  //get location and transform it
269  double3 &gip_pt = sub_gip.gip_points[gip_inx];
270 
271  //get value of a shape function
272  Scalar shape_value0 = sub_shape.svals[H2D_HCFE_VALUE0][gip_inx];
273  Scalar shape_value1 = sub_shape.svals[H2D_HCFE_VALUE1][gip_inx];
274  Scalar shape_curl = sub_shape.svals[H2D_HCFE_CURL][gip_inx];
275 
276  //get value of ref. solution
277  Scalar ref_value0 = sub_trf.coef_mx * rval[son][H2D_HCFE_VALUE0][gip_inx];
278  Scalar ref_value1 = sub_trf.coef_my * rval[son][H2D_HCFE_VALUE1][gip_inx];
279  //coef_curl * curl
280  Scalar ref_curl = coef_curl * rval[son][H2D_HCFE_CURL][gip_inx];
281 
282  //evaluate a right-hand value
283  Scalar value = (shape_value0 * ref_value0)
284  + (shape_value1 * ref_value1)
285  + (shape_curl * ref_curl);
286 
287  total_value += gip_pt[H2D_GIP2D_W] * value;
288  }
289  return total_value;
290  }
291 
292  template<typename Scalar>
293  double HcurlProjBasedSelector<Scalar>::evaluate_error_squared_subdomain(Element* sub_elem, const typename ProjBasedSelector<Scalar>::ElemGIP& sub_gip, int son, const typename ProjBasedSelector<Scalar>::ElemSubTrf& sub_trf, const typename ProjBasedSelector<Scalar>::ElemProj& elem_proj, Scalar* rval[H2D_MAX_ELEMENT_SONS][MAX_NUMBER_FUNCTION_VALUES_FOR_SELECTORS])
294  {
295  double total_error_squared = 0;
296  double coef_curl = std::abs(sub_trf.coef_mx * sub_trf.coef_my);
297  for (int gip_inx = 0; gip_inx < sub_gip.num_gip_points; gip_inx++)
298  {
299  //get location and transform it
300  double3 &gip_pt = sub_gip.gip_points[gip_inx];
301 
302  //calculate value of projected solution
303  Scalar proj_value0 = 0, proj_value1 = 0, proj_curl = 0;
304  for (int i = 0; i < elem_proj.num_shapes; i++)
305  {
306  int shape_inx = elem_proj.shape_inxs[i];
307  proj_value0 += elem_proj.shape_coeffs[i] * elem_proj.svals[shape_inx][H2D_HCFE_VALUE0][gip_inx];
308  proj_value1 += elem_proj.shape_coeffs[i] * elem_proj.svals[shape_inx][H2D_HCFE_VALUE1][gip_inx];
309  proj_curl += elem_proj.shape_coeffs[i] * elem_proj.svals[shape_inx][H2D_HCFE_CURL][gip_inx];
310  }
311 
312  {
313  //get value of ref. solution
314  Scalar ref_value0 = sub_trf.coef_mx * rval[son][H2D_HCFE_VALUE0][gip_inx];
315  Scalar ref_value1 = sub_trf.coef_my * rval[son][H2D_HCFE_VALUE1][gip_inx];
316  //coef_curl * curl
317  Scalar ref_curl = coef_curl * rval[son][H2D_HCFE_CURL][gip_inx];
318 
319  //evaluate error
320  double error_squared = sqr(proj_value0 - ref_value0)
321  + sqr(proj_value1 - ref_value1)
322  + sqr(proj_curl - ref_curl);
323 
324  total_error_squared += gip_pt[H2D_GIP2D_W] * error_squared;
325  }
326  }
327  return total_error_squared;
328  }
329  template class HERMES_API HcurlProjBasedSelector < double > ;
330  template class HERMES_API HcurlProjBasedSelector < std::complex<double> > ;
331  }
332  }
333 }
double coef_my
A coefficient that scales df/dy for each subdomain. A coefficient represents effects of a transformat...
Element * get_active_element() const
Definition: transformable.h:48
virtual void precalc_shapes(const double3 *gip_points, const int num_gip_points, const Trf *trfs, const int num_noni_trfs, const std::vector< typename OptimumSelector< Scalar >::ShapeInx > &shapes, const int max_shape_inx, typename ProjBasedSelector< Scalar >::TrfShape &svals, ElementMode2D mode)
Calculates values of shape function at GIP for all transformations.
Definition: adapt.h:24
virtual void precalc_ref_solution(int inx_son, MeshFunction< Scalar > *rsln, Element *element, int intr_gip_order, Scalar *rval[H2D_MAX_ELEMENT_SONS][MAX_NUMBER_FUNCTION_VALUES_FOR_SELECTORS])
Returns an array of values of the reference solution at integration points.
A projection-based selector for Hcurl space.
Definition: function.h:36
Stores one element of a mesh.
Definition: element.h:107
virtual void get_current_order_range(Element *element, int &min_order, int &max_order)
Sets OptimumSelector::current_max_order and OptimumSelector::current_min_order.
virtual void free_ref_solution_data(int inx_son, Scalar *rval[H2D_MAX_ELEMENT_SONS][MAX_NUMBER_FUNCTION_VALUES_FOR_SELECTORS])
Frees the data allocated in precalc_ref_solution.
Represents a function defined on a mesh.
Definition: mesh_function.h:56
Quad2D * get_quad_2d() const
Returns the current quadrature points.
Definition: function.cpp:129
Integration points in the reference domain of an element of a candidate.
double coef_mx
A coefficient that scales df/dx for each subdomain. A coefficient represents effects of a transformat...
Index of df/dx.
Definition: shapeset.h:27
virtual const Scalar * get_dx_values(int component=0) const
Returns the x partial derivative.
Definition: function.cpp:175
virtual const Scalar * get_dy_values(int component=0) const
Returns the y partial derivative.
Definition: function.cpp:182
Index of df/dy.
Definition: shapeset.h:28
double3 * gip_points
Integration points and weights. The first index is an index of an integration point, the second index is defined through the enum GIP2DIndices.
A general projection-based selector.
Definition: function.h:33
CandList
Predefined list of candidates.
Definition: candidates.h:46
virtual double evaluate_error_squared_subdomain(Element *sub_elem, const typename ProjBasedSelector< Scalar >::ElemGIP &sub_gip, int son, const typename ProjBasedSelector< Scalar >::ElemSubTrf &sub_trf, const typename ProjBasedSelector< Scalar >::ElemProj &elem_proj, Scalar *rval[H2D_MAX_ELEMENT_SONS][MAX_NUMBER_FUNCTION_VALUES_FOR_SELECTORS])
Evaluates an squared error of a projection of an element of a candidate onto subdomains.
Index of a function value f.
Definition: shapeset.h:26
unsigned short iro_cache
Increase in integration order, see RefMap::calc_inv_ref_order()
Definition: element.h:199
Y-axis coordinate.
Definition: quad.h:29
Scalar * shape_coeffs
Coefficients of shape indices of a projection.
double2 t
The 2x2 diagonal transformation matrix.
Definition: transformable.h:32
2D transformation.
Definition: transformable.h:29
H(curl) shapeset with Legendre bubbles and gradients of H1 functions as edges.
std::vector< TrfShapeExp > & svals
A precalculated shape-function values. Empty is not defined.
virtual Scalar evaluate_rhs_subdomain(Element *sub_elem, const typename ProjBasedSelector< Scalar >::ElemGIP &sub_gip, int son, const typename ProjBasedSelector< Scalar >::ElemSubTrf &sub_trf, const typename ProjBasedSelector< Scalar >::ElemSubShapeFunc &sub_shape, Scalar *rval[H2D_MAX_ELEMENT_SONS][MAX_NUMBER_FUNCTION_VALUES_FOR_SELECTORS])
Evaluates a value of the right-hande side in a subdomain.
TrfShapeExp & svals
Evaluate values of a shape function. If TrfShapeExp::empty(), no precalculated values are available...
A transformation from a reference domain of a subdomain to a reference domain of an element of a cand...
virtual void precalc_ortho_shapes(const double3 *gip_points, const int num_gip_points, const Trf *trfs, const int num_noni_trfs, const std::vector< typename OptimumSelector< Scalar >::ShapeInx > &shapes, const int max_shape_inx, typename ProjBasedSelector< Scalar >::TrfShape &svals, ElementMode2D mode)
Calculates values of orthogonalized shape function at GIP for all transformations.
#define H2D_MAX_ELEMENT_SONS
Macros.
Definition: global.h:30
virtual const Scalar * get_fn_values(int component=0) const
Returns function values.
Definition: function.cpp:168
virtual double ** build_projection_matrix(double3 *gip_points, int num_gip_points, const int *shape_inx, const int num_shapes, ElementMode2D mode)
Builds projection matrix using a given set of shapes.
HcurlProjBasedSelector(CandList cand_list=H2D_HP_ANISO, int max_order=H2DRS_DEFAULT_ORDER, HcurlShapeset *user_shapeset=nullptr)
Constructor.
#define H2DRS_DEFAULT_ORDER
A default order. Used to indicate an unkonwn order or a maximum support order.
Definition: global.h:103
X-axis coordinate.
Definition: quad.h:27