Hermes2D  3.0
l2_proj_based_selector.cpp
1 #include "matrix.h"
2 #include "l2_proj_based_selector.h"
3 
4 namespace Hermes
5 {
6  namespace Hermes2D
7  {
8  namespace RefinementSelectors
9  {
10  template<typename Scalar>
11  L2ProjBasedSelector<Scalar>::L2ProjBasedSelector(CandList cand_list, int max_order, L2Shapeset* user_shapeset)
12  : ProjBasedSelector<Scalar>(cand_list, max_order, user_shapeset == nullptr ? new L2Shapeset() : user_shapeset, Range(1, 1), Range(0, H2DRS_MAX_ORDER)), user_shapeset(user_shapeset == nullptr ? false : true)
13  {
14  if (user_shapeset != nullptr)
15  {
16  this->warn("Warning: The user shapeset provided for the selector has to have a correct copy constructor implemented.");
17  this->warn("Warning: The functionality for cloning user shapeset is to be implemented yet.");
18  }
19  }
20 
21  template<typename Scalar>
23  {
24  if (!this->user_shapeset)
25  delete this->shapeset;
26  }
27 
28  template<typename Scalar>
29  void L2ProjBasedSelector<Scalar>::get_current_order_range(Element* element, int& min_order_, int& max_order_)
30  {
31  int max_element_order = (20 - element->iro_cache) / 2 - 1;
32  if (this->max_order == H2DRS_DEFAULT_ORDER)
33  max_order_ = max_element_order;
34  else
35  max_order_ = std::min(this->max_order, max_element_order);
36  min_order_ = 0;
37  }
38 
39  template<typename Scalar>
40  void L2ProjBasedSelector<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)
41  {
42  //for all transformations
43  bool done = false;
44  int inx_trf = 0;
45  while (!done && inx_trf < H2D_TRF_NUM)
46  {
47  //prepare data for processing
48  const Trf& trf = trfs[inx_trf];
49  std::vector<typename ProjBasedSelector<Scalar>::TrfShapeExp>& trf_svals = svals[inx_trf];
50 
51  //allocate
52  trf_svals.resize(max_shape_inx + 1);
53 
54  //for all shapes
55  const int num_shapes = (int)shapes.size();
56  for (int i = 0; i < num_shapes; i++)
57  {
58  int inx_shape = shapes[i].inx;
59  typename ProjBasedSelector<Scalar>::TrfShapeExp& shape_exp = trf_svals[inx_shape];
60 
61  //allocate
62  shape_exp.allocate(H2D_L2FE_NUM, num_gip_points);
63 
64  //for all GIP points
65  for (int k = 0; k < num_gip_points; k++)
66  {
67  //transform coordinates
68  double ref_x = gip_points[k][H2D_GIP2D_X] * trf.m[0] + trf.t[0];
69  double ref_y = gip_points[k][H2D_GIP2D_Y] * trf.m[1] + trf.t[1];
70 
71  //for all expansions: retrieve values
72  shape_exp[H2D_L2FE_VALUE][k] = this->shapeset->get_fn_value(inx_shape, ref_x, ref_y, 0, mode);
73  }
74  }
75 
76  //move to the next transformation
77  if (inx_trf == H2D_TRF_IDENTITY)
78  done = true;
79  else
80  {
81  inx_trf++;
82  if (inx_trf >= num_noni_trfs) //if all transformations were processed, move to the identity transformation
83  inx_trf = H2D_TRF_IDENTITY;
84  }
85  }
86  if (!done)
87  throw Exceptions::Exception("All transformation processed but identity transformation not found.");
88  }
89 
90  template<typename Scalar>
91  void L2ProjBasedSelector<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)
92  {
93  //calculate values
94  precalc_shapes(gip_points, num_gip_points, trfs, num_noni_trfs, shapes, max_shape_inx, svals, mode);
95 
96  //calculate orthonormal basis
97  const int num_shapes = (int)shapes.size();
98  for (int i = 0; i < num_shapes; i++)
99  {
100  const int inx_shape_i = shapes[i].inx;
101 
102  //orthogonalize
103  for (int j = 0; j < i; j++)
104  {
105  const int inx_shape_j = shapes[j].inx;
106 
107  //calculate product of non-transformed functions
108  double product = 0.0;
109  for (int k = 0; k < num_gip_points; k++)
110  {
111  double sum = 0.0;
112  sum += svals[H2D_TRF_IDENTITY][inx_shape_i][H2D_L2FE_VALUE][k] * svals[H2D_TRF_IDENTITY][inx_shape_j][H2D_L2FE_VALUE][k];
113  product += gip_points[k][H2D_GIP2D_W] * sum;
114  }
115 
116  //for all transformations
117  int inx_trf = 0;
118  bool done = false;
119  while (!done && inx_trf < H2D_TRF_NUM)
120  {
121  //for all integration points
122  for (int k = 0; k < num_gip_points; k++)
123  {
124  svals[inx_trf][inx_shape_i][H2D_L2FE_VALUE][k] -= product * svals[inx_trf][inx_shape_j][H2D_L2FE_VALUE][k];
125  }
126 
127  //move to the next transformation
128  if (inx_trf == H2D_TRF_IDENTITY)
129  done = true;
130  else
131  {
132  inx_trf++;
133  if (inx_trf >= num_noni_trfs) //if all transformations were processed, move to the identity transformation
134  inx_trf = H2D_TRF_IDENTITY;
135  }
136  }
137  //identity transformation has to be the last transformation
138  if (!done)
139  throw Exceptions::Exception("All transformation processed but identity transformation not found.");
140  }
141 
142  //normalize
143  //calculate norm
144  double norm_squared = 0.0;
145  for (int k = 0; k < num_gip_points; k++)
146  {
147  double sum = 0.0;
148  sum += sqr(svals[H2D_TRF_IDENTITY][inx_shape_i][H2D_L2FE_VALUE][k]);
149  norm_squared += gip_points[k][H2D_GIP2D_W] * sum;
150  }
151  double norm = sqrt(norm_squared);
152  if (!finite(1 / norm))
153  throw Exceptions::Exception("Norm (%g) is almost zero.", norm);
154 
155  //for all transformations: normalize
156  int inx_trf = 0;
157  bool done = false;
158  while (!done && inx_trf < H2D_TRF_NUM)
159  {
160  //for all integration points
161  for (int k = 0; k < num_gip_points; k++)
162  {
163  svals[inx_trf][inx_shape_i][H2D_L2FE_VALUE][k] /= norm;
164  }
165 
166  //move to the next transformation
167  if (inx_trf == H2D_TRF_IDENTITY)
168  done = true;
169  else
170  {
171  inx_trf++;
172  if (inx_trf >= num_noni_trfs) //if all transformations were processed, move to the identity transformation
173  inx_trf = H2D_TRF_IDENTITY;
174  }
175  }
176  //identity transformation has to be the last transformation
177  if (!done)
178  throw Exceptions::Exception("All transformation processed but identity transformation not found.");
179  }
180  }
181 
182  template<typename Scalar>
183  void L2ProjBasedSelector<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])
184  {
185  const int num_gip = rsln->get_quad_2d()->get_num_points(intr_gip_order, rsln->get_active_element()->get_mode());
186  rval[inx_son][H2D_L2FE_VALUE] = malloc_with_check<Scalar>(num_gip);
187  memcpy(rval[inx_son][H2D_L2FE_VALUE], rsln->get_fn_values(0), num_gip * sizeof(Scalar));
188  }
189 
190  template<typename Scalar>
191  void L2ProjBasedSelector<Scalar>::free_ref_solution_data(int inx_son, Scalar* rval[H2D_MAX_ELEMENT_SONS][MAX_NUMBER_FUNCTION_VALUES_FOR_SELECTORS])
192  {
193  free_with_check(rval[inx_son][H2D_L2FE_VALUE]);
194  }
195 
196  template<typename Scalar>
197  double** L2ProjBasedSelector<Scalar>::build_projection_matrix(double3* gip_points, int num_gip_points,
198  const int* shape_inx, const int num_shapes, ElementMode2D mode)
199  {
200  //allocate
201  double** matrix = new_matrix<double>(num_shapes, num_shapes);
202 
203  //calculate products
204  int inx_row = 0;
205  for (int i = 0; i < num_shapes; i++, inx_row += num_shapes)
206  {
207  double* matrix_row = matrix[i];
208  int shape0_inx = shape_inx[i];
209  for (int k = 0; k < num_shapes; k++)
210  {
211  int shape1_inx = shape_inx[k];
212 
213  double value = 0.0;
214  for (int j = 0; j < num_gip_points; j++)
215  {
216  double gip_x = gip_points[j][H2D_GIP2D_X], gip_y = gip_points[j][H2D_GIP2D_Y];
217  double value0 = this->shapeset->get_value(H2D_FEI_VALUE, shape0_inx, gip_x, gip_y, 0, mode);
218  double value1 = this->shapeset->get_value(H2D_FEI_VALUE, shape1_inx, gip_x, gip_y, 0, mode);
219 
220  value += gip_points[j][H2D_GIP2D_W] * (value0*value1);
221  }
222 
223  matrix_row[k] = value;
224  }
225  }
226 
227  return matrix;
228  }
229 
230  template<typename Scalar>
231  Scalar L2ProjBasedSelector<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])
232  {
233  Scalar total_value = 0;
234  for (int gip_inx = 0; gip_inx < sub_gip.num_gip_points; gip_inx++)
235  {
236  //get location and transform it
237  double3 &gip_pt = sub_gip.gip_points[gip_inx];
238 
239  //get value of a shape function
240  double shape_value = sub_shape.svals[H2D_L2FE_VALUE][gip_inx];
241 
242  //get value of ref. solution
243  Scalar ref_value;
244  ref_value = rval[son][H2D_L2FE_VALUE][gip_inx];
245 
246  //evaluate a right-hand value
247  Scalar value = (shape_value * ref_value);
248 
249  total_value += gip_pt[H2D_GIP2D_W] * value;
250  }
251  return total_value;
252  }
253 
254  template<typename Scalar>
255  double L2ProjBasedSelector<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])
256  {
257  double total_error_squared = 0;
258  for (int gip_inx = 0; gip_inx < sub_gip.num_gip_points; gip_inx++)
259  {
260  //get location and transform it
261  double3 &gip_pt = sub_gip.gip_points[gip_inx];
262 
263  //calculate value of projected solution
264  Scalar proj_value = 0;
265  for (int i = 0; i < elem_proj.num_shapes; i++)
266  {
267  int shape_inx = elem_proj.shape_inxs[i];
268  proj_value += elem_proj.shape_coeffs[i] * elem_proj.svals[shape_inx][H2D_L2FE_VALUE][gip_inx];
269  }
270 
271  //get value of ref. solution
272  Scalar ref_value = rval[son][H2D_L2FE_VALUE][gip_inx];
273 
274  //evaluate error
275  double error_squared = sqr(proj_value - ref_value);
276 
277  total_error_squared += gip_pt[H2D_GIP2D_W] * error_squared;
278  }
279 
280  return total_error_squared;
281  }
282  template class HERMES_API L2ProjBasedSelector < double > ;
283  template class HERMES_API L2ProjBasedSelector < std::complex<double> > ;
284  }
285  }
286 }
virtual void get_current_order_range(Element *element, int &min_order, int &max_order)
Sets OptimumSelector::current_max_order and OptimumSelector::current_min_order.
Element * get_active_element() const
Definition: transformable.h:48
Definition: adapt.h:24
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.
Stores one element of a mesh.
Definition: element.h:107
#define H2DRS_MAX_ORDER
A maximum order suported by refinement selectors.
Definition: global.h:104
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.
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.
L2 shapeset - products of legendre polynomials.
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
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
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.
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
std::vector< TrfShapeExp > & svals
A precalculated shape-function values. Empty is not defined.
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...
L2ProjBasedSelector(CandList cand_list=H2D_HP_ANISO, int max_order=H2DRS_DEFAULT_ORDER, L2Shapeset *user_shapeset=nullptr)
Constructor.
#define H2D_MAX_ELEMENT_SONS
Macros.
Definition: global.h:30
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.
virtual const Scalar * get_fn_values(int component=0) const
Returns function values.
Definition: function.cpp:168
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.
A projection-based selector for L2 space.
Definition: function.h:34
#define H2DRS_DEFAULT_ORDER
A default order. Used to indicate an unkonwn order or a maximum support order.
Definition: global.h:103
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.
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.
X-axis coordinate.
Definition: quad.h:27