refinement.hh

Go to the documentation of this file.
1 #ifndef refinement_hh
2 #define refinement_hh
3 
4 //#include "estimator.hh"
5 #include "estimator/marking.hh"
6 #include "space/hpMethod.hh"
8 #include "hp2D/quad.hh"
10 //#include "lInfinityNorm.hh"
11 //#include "geometry/diameter.hh"
12 
13 #if __GNUC__ == 2
14 #include <float.h>
15 #define EPS DBL_EPSILON
16 #define MAX_D DBL_MAX
17 #else
18 #include <limits>
19 #define MAX_D std::numeric_limits<double>::max()
20 #endif
21 
22 using concepts::Real;
23 
24 namespace estimator{
25 
29 template<class F=Real>
31 
32 public:
33  //apply routine that adjust the elements in the space,
34  //e.g. elements that are marked for refinement
35 
36  virtual void apply(concepts::SpaceOnCells<F>& space) = 0;
37 
38  protected:
39  virtual std::ostream& info(std::ostream& os) const = 0;
40 
41 };
42 
43 
48 //template<class F=Real>
49 //class H_Refinement: public concepts::Refinement<F>{
50 //public:
51 // H_Refinement(const concepts::Marking& mark): mark_(mark){}
52 //
53 //private :
54 // concepts::Marking& mark_;
55 //};
56 
57 // Abstract Class for hp refinement strategies.
58 template<class F = Real, ushort dim = 2>
59 class HP_Refinement : public Refinement<F>{
60 public :
61 
62  virtual ~HP_Refinement(){};
63 
64  virtual void apply(concepts::SpaceOnCells<F>& space){
65  hp2D::hpAdaptiveSpace<F>* spc = dynamic_cast<hp2D::hpAdaptiveSpace<F>*>(&space);
66  if(spc){
67  std::unique_ptr<hp2D::hpAdaptiveSpaceH1::Scan> sc(spc->scan());
68  while (!sc->eos()) {
69  concepts::ElementWithCell<Real>& element = (*sc)++;
70  const hp2D::Quad<Real>* quad = dynamic_cast<const hp2D::Quad<Real>*> (&element);
71  uint K = quad->support().key();
72 
73  typename concepts::HashMap<concepts::AdaptiveAdjustP<dim> >::const_iterator iter;
74  //current element has a refinement rule
75  if((iter = adjusts_.find(K)) != adjusts_.end())
76  spc->adjust(*quad, iter->second);
77  }
78  spc->rebuild();
79  }else{
80  throw conceptsException(concepts::MissingFeature("Currently just hpAdaptiveSpace supported"));
81  }
82 
83  }
84 
85 protected :
86  //map from cellKey to its adjust strategy
88 
89 };
90 
98 template<class F = Real, ushort dim = 2>
99 class Prediction : public HP_Refinement<F,dim>{
100 public:
101 
102  //starting prediction assumption
103  // if cell was not refined before :
104  // H : h refinement is forced for first refinement
105  // P : p refinement is forced for first refinement
106  // after prediction computation is possible
107  enum Mode {H = 0 , P = 1};
108 
122  Prediction(const concepts::SpaceOnCells<F>& space, enum Mode mode = P,
123  const Real hdec = sqrt(4.0), const Real pdec = sqrt(0.4), const Real nRdec = 1.0):
124  hdecay_(hdec), pdecay_(pdec), noRefdecay_(nRdec)
125  {
126  const hp2D::hpAdaptiveSpace<F>* spc = dynamic_cast<const hp2D::hpAdaptiveSpace<F>* >(&space);
127  if(spc){
128  //information about the polynomial degree
129  hp2D::hpFull& hpfull = spc->prebuild();
130  std::unique_ptr<hp2D::hpAdaptiveSpaceH1::Scan> sc(spc->scan());
131  while (*sc) {
132  concepts::ElementWithCell<Real>& elm = (*sc)++;
133  //key of current element
134  uint K = elm.cell().connector().key().key();
135  //set default prediction strategy according to chosen mode
136  pred_[K] = mode ? MAX_D : 0;
137  //add cell information
138  cells_[K] = &elm.cell();
139 
140  const concepts::Connector2* cntr =
141  dynamic_cast<const concepts::Connector2*>(&elm.cell().connector());
142  //the cntr must be 2d since it is out of hpAdaptiveSpace
144  //set maximal polynomial degree (i.e inner), if no inner exists
145  // default is 1 if there is a vertex dof on this cell
146  const ushort* pMax = hpfull.innerDof(*cntr);
147  //no inner dof
148  if(pMax == 0){
149  //check if there is at least one dof, i.e a vertex dof
150  //this has 4 vertexList searches (costs?)
151  bool hasVDof = (hpfull.vtxDof(*cntr, 0) || hpfull.vtxDof(*cntr, 1) ||
152  hpfull.vtxDof(*cntr, 2) || hpfull.vtxDof(*cntr, 3));
153  if(hasVDof){
154  //set default to one
155  p_[K] = concepts::Array<ushort>(dim);
156  for(ushort i = 0 ; i < dim; ++i)
157  p_[K][i] = 1;
158  }else{
159  //no dofs on vertices, e.g. :
160  //
161  // Cell x for h=1 only has hanging nodes as vertices and is therefore not supported
162  // __________
163  // |__________|
164  // | |x| |
165  // | |_|
166  // |___|_|____|
167  throw conceptsException(concepts::MissingFeature("Cell has no vertex dofs, cannot be handled atm"));
168  }
169  }else{
170  //has inner Dof and therefore a polynomial degree >= 2
171  p_[K] = concepts::Array<ushort>(dim);
172  for(ushort i = 0 ; i < dim; ++i)
173  p_[K][i] = pMax[i];
174  }
175 
176  }
177  }else{
178  throw conceptsException(concepts::MissingFeature("Currently just hpAdaptiveSpace supported"));
179  }
180  }
181 
182 
183  virtual ~Prediction(){}
184 
185 
186 // builds the refinement strategys for given mark and estimator
188 
189  //all marked Cells (its keys)
190  const concepts::Set<uint>& marked = mark.marks();
191 
192 
193  //check all cells if they are marked
195  for( ; cIter != cells_.end(); ++cIter){
196  //current cell
197  uint K = cIter->first;
198  //cell is marked for refinement
199  if(marked.exist(K)){
200  //DEBUGL(1, "on K = "<< K << " error = "<< estimator(K) << " pred = "<< pred_[K]);
201  if (estimator(K) > pred_[K]){
202  //DEBUGL(1, "on K = "<< K << " error = "<< estimator(K) << " pred = "<< pred_[K] << " apply h ref");
203  this->adjusts_[K] = concepts::AdaptiveAdjustP<dim>(1, 0);
204  //add cell to as href marked, later needed to compute child prediction
205  //those are not avaiable atm
206  hRefCurCells_.push_back(cIter->second);
207  //only isotropic p refinement is allowed
208  conceptsAssert(p_[K][0]==p_[K][1], concepts::Assertion());
209  //cell will have childs later, now we save child info for it here
210  pred_[K] = 0.5 * hdecay_ * std::pow(0.5,p_[K][0] ) * estimator(K);
211  //also polynomial degree is setted on kids etc
212 
213  //DEBUGL(1, "set prediction for error = "<< estimator(K) << " mult = "<< 0.5 * hdecay_ * std::pow(0.5,p_[K][0] )
214  // << " so pred = " << pred_[K]);
215  }
216  //set for p refinement
217  else{
218  //DEBUGL(1, "on K = "<< K << " error = "<< estimator(K) << " pred = "<< pred_[K] << " apply p ref");
219  //add isotropic p refinement plus 1
220  this->adjusts_[K] = concepts::AdaptiveAdjustP<dim>(0, 1);
221  //update info of polynomial degree on this element
222  p_[K][0] += 1; p_[K][1] += 1;
223  //DEBUGL(1, "old pred_[K] = "<< pred_[K] << " pdecay_ = "<< pdecay_ << " new :");
224  pred_[K] = pdecay_ * estimator(K);
225  }
226 
227  //DEBUGL(1, "prediction( K = "<< K << " ) "<< pred_[K]);
228 
229  }else{
230  //not marked update with no refinement decay
231  pred_[K] = noRefdecay_ * pred_[K];
232  }
233  }
234 }
235 
236 
237 virtual void apply(concepts::SpaceOnCells<F>& space){
238  //refine space
240  // clear all adjusts as they needed to be evaluated new for same elements
241  this->adjusts_.erase(this->adjusts_.begin(), this->adjusts_.end());
242 
243  //update informations for predictions, i.e. for h refined elements, that after apply
244  //have new children
246  for(cIter = hRefCurCells_.begin(); cIter != hRefCurCells_.end(); ++cIter){
247 
248  //current adjacent cell key
249  uint K = (*cIter)->connector().key().key();
250  const concepts::Cell* child = 0;
251  //current child it must be 4, else not supported
252  //since prediction must be adapted for halfsplitting, e.g. in 2D
253  uint i = 0 ;
254  for( ; (child = (*cIter)->child(i)) !=0 ; ++i){
255  //key of builded child cell
256  uint cK = child->connector().key();
257  //update polynomial degree to child cells
258  p_[cK] = concepts::Array<ushort>(dim);
259  for(ushort i = 0 ; i < dim; ++i)
260  p_[cK][i] = p_[K][i];
261 
262  pred_[cK] = pred_[K];
263  //add new cell to cell list
264  cells_[cK] = child;
265  }
266  //erase old polynomial degree information
267  p_.erase(p_.find(K));
268  //erase old prediction info
269  pred_.erase(pred_.find(K));
270  //erase coarse cell
271  cells_.erase(cells_.find(K));
272 
273 
274  }
275  //clear this map, since those cells no longer are finest cells in space
276  hRefCurCells_.erase(hRefCurCells_.begin(), hRefCurCells_.end());
277 
278 
279 }
280 
281 protected:
282  virtual std::ostream& info(std::ostream& os) const{
283  os << " Prediction[ type = hprefinement Strategy ]" << std::endl;
284  os << "refinement chooses : "<< std::endl;
285  return os << this->adjusts_;
286  //further output : history
287  // the parameters etc
288  }
289 
290 
291 private:
292  //error decay value for h refinement
294  //error decay value for p refinement
296  //error decay value for when no is applied refinement
298 
299  //map from cell Key to prediction
301 
302  //map from cellKey to its highest polynomial degree, atm just isotropic, therefore only one ushort
304 
305  //current cells of given space
307 
308  //map from (static) key to belonging that is hrefined in the space and therefore has kid elements
309  //that need updated prediction values in the moment when the space is rebuilded
310  //current and updated with build / apply routine
312 
313 
314 };
315 
316 
317 
318 
319 
320 
321 
322 //template<class F = Real, ushort dim = 2>
323 //class Embedding : public HP_Refinement<F,dim>{
324 //public:
325 //
326 // // H1 just H1 norm
327 // // H2Partial2 is H2 norm with just dxx and dyy part of H2 seminorm
328 // enum indicator {H1=0 , H2Partial2};
329 //
331 // *
332 // * @param reg regularity barrier parameter that decides which refinement should be applied
333 // * i.e for indicator smaller as reg we apply h-refinement
334 // * else p refinement
335 // */
336 // Embedding(const concepts::SpaceOnCells<F>& space, const concepts::Vector<F>& sol,const Real reg, enum indicator ind = H2Partial2)
337 // :spc_(space), sol_(sol), reg_(reg), ind_(ind)
338 // {
339 //
340 // const hp2D::hpAdaptiveSpace<F>* spc = dynamic_cast<const hp2D::hpAdaptiveSpace<F>* >(&space);
341 // if(spc){
342 // //information about the polynomial degree
343 // hp2D::hpFull& hpfull = spc->prebuild();
344 // std::auto_ptr<hp2D::hpAdaptiveSpaceH1::Scan> sc(spc->scan());
345 // while (*sc) {
346 // concepts::ElementWithCell<Real>& elm = (*sc)++;
347 //
348 // const hp2D::Quad<Real>* quad = dynamic_cast<const hp2D::Quad<Real>*>(&elm);
349 // conceptsAssert(quad, concepts::Assertion());
350 // //key of current element
351 // uint K = elm.cell().connector().key().key();
352 // //add cell information
353 // quads_[K] = quad;
354 // }
355 // }else{
356 // throw conceptsException(concepts::MissingFeature("Currently just hpAdaptiveSpace supported"));
357 // }
358 // }
359 //
360 //
361 //
362 //void buildRefinement(const Marking<F>& mark, const concepts::LocalEstimator<F>& estimator){
363 //
364 //
365 // //all marked Cells (its keys)
366 // const concepts::Set<uint>& marked = mark.marks();
367 //
368 // //diameter map from each Element /TODO: Just take part via Set constructor
369 // concepts::Diameter h_K(spc_,concepts::Diameter::PLANE,true);
370 //
371 // //potential used function evaluations
372 // concepts::ElementFormulaVector<1> apprxSol(spc_, sol_, hp2D::Value<Real>());
373 // concepts::ElementFormulaVector<2> apprxGrad(spc_, sol_, hp2D::Grad<Real>());
374 // concepts::ElementFormulaVector<1> partial_xx(spc_, sol_, hp2D::Partial_xx<Real,hp2D::Quad<Real> >());
375 // concepts::ElementFormulaVector<1> partial_yy(spc_, sol_, hp2D::Partial_yy<Real,hp2D::Quad<Real> >());
376 //
377 // //check all cells if they are marked
378 // concepts::Set<uint>::const_iterator mIter = marked.begin();
379 // for( ; mIter != marked.end(); ++mIter){
380 // //current marked cell key
381 // uint K = *mIter;
382 // //cell is marked for refinement
383 //
384 // const hp2D::Quad<Real>* quad = quads_[K];
385 //
386 // Real indicator = 0;
387 // //TODO: build static change of point rate and method
388 // estimator::LInfinity maximum(*quad, sol_,1,2);
389 // // L^2(K) norm squared
390 // Real L_2_K = concepts::L2product(*quad, apprxSol);
391 // // H1 seminorm squared
392 // Real H_1_K = concepts::L2product(*quad, apprxGrad);
393 //
394 //
395 // if(ind_ == H1){
396 // indicator = std::pow(maximum.get(),2) / ( std::pow( 0.5 * h_K(K), -2) * L_2_K + H_1_K );
397 // //DEBUGL(1, "ind H (K = " << K <<") = "<< indicator);
398 // }
399 // else if (ind_ == H2Partial2){
400 // // H2 partial xx seminorm
401 // Real H_2xx = concepts::L2product(*quad, partial_xx);
402 // // H2 partial yy seminorm
403 // Real H_2yy = concepts::L2product(*quad, partial_yy);
404 // indicator = std::pow(maximum.get(), 2) / (std::pow(
405 // 0.5 * h_K(K), -2) * L_2_K + H_1_K + std::pow(0.5 * h_K(
406 // K), 2) * (H_2xx + H_2yy));
407 // //DEBUGL(1, "ind G (K = " << K <<") = "<< indicator);
408 // }else{
409 // throw conceptsException(concepts::MissingFeature("Other indicator types not implemented"));
410 // }
411 //
412 //
413 // // set adaptive strategy
414 // if( indicator < reg_ ){
415 // this->adjusts_[K] = concepts::AdaptiveAdjustP<dim>(1, 0);
416 // DEBUGL(1, "ind (K = " << K <<") = "<< indicator << " --> h-refine");
417 // }
418 // else{
419 // this->adjusts_[K] = concepts::AdaptiveAdjustP<dim>(0, 1);
420 // DEBUGL(1, "ind (K = " << K <<") = "<< indicator << " --> p-refine");
421 // }
422 // }//loop over all marked cells
423 //}//build refinement
424 //
425 //protected:
426 // virtual std::ostream& info(std::ostream& os) const{
427 // os << " Embedding[ type = hp - Refinement Strategy :" << std::endl;
428 // os << "Indicator type :"<< ( (ind_) ? "H2Partial2" : "H1" ) << std::endl;
429 // return os << "regularity barrier = "<< reg_ << "]";
430 // //further output : history
431 // // the parameters etc
432 // }
433 //
434 //
435 //
436 //private:
437 // const concepts::SpaceOnCells<F>& spc_;
438 //
439 // const concepts::Vector<F>& sol_;
440 //
441 // concepts::HashMap<const hp2D::Quad<Real>* > quads_;
442 // //regularity barrier
443 // Real reg_;
444 //
445 // //indicator type
446 // enum indicator ind_;
447 //
448 //};
449 
450 
451 
452 
453 
454 
455 
456 
457 
458 
459 
465 template<class F = Real, ushort dim = 2>
466 class AprioriVertex : public HP_Refinement<F,dim>{
467 public:
468 
473  AprioriVertex(const concepts::SpaceOnCells<F>& space, const uint attrb){
474  attrbs_.insert(attrb);
475  //this easily extends to non hp2D::hpAdaptiveSpace
476  std::unique_ptr<typename concepts::SpaceOnCells<F>::Scanner> sc(space.scan());
477  while (*sc) {
478  concepts::ElementWithCell<F>& elm = (*sc)++;
479  //key of current element
480  uint K = elm.cell().connector().key().key();
481  //add cell information
482  cells_[K] = &elm.cell();
483  }
484  }
485 
486 
487  virtual ~AprioriVertex(){}
488 
489  void addVertex(const uint attrb){
490  attrbs_.insert(attrb);
491  }
492 
493  //to generalize, here estimator is not needed
495  const concepts::Set<uint>& marks = mark.marks();
496  for(concepts::Set<uint>::const_iterator mIter = marks.begin(); mIter != marks.end(); ++mIter ){
497  //key of marked element
498  uint K = *mIter;
499 
500 // this->adjusts_[K] = concepts::AdaptiveAdjustP<dim>(1, 0);
501 
502  const concepts::Connector2* cntr =
503  dynamic_cast<const concepts::Connector2*>(&(cells_[K]->connector()));
504  if(cntr){
505  bool pRef = true;
506  //current vertex
507  const concepts::Connector0* vertex = 0;
508  for( uint i = 0; (vertex = cntr->vertex(i) ) != 0; ++i){
509  //quad has a singular marked vertex
510  if(attrbs_.exist(vertex->attrib().attrib())){
511  this->adjusts_[K] = concepts::AdaptiveAdjustP<dim>(1, 0);
512  pRef = false;
513  break;
514  }
515  }
516  //if no h refinement is applied
517  if(pRef)
518  this->adjusts_[K] = concepts::AdaptiveAdjustP<dim>(0, 1);
519 
520  }else{
521  throw conceptsException(concepts::MissingFeature("Currently Only 2d elements supported"));
522  }
523  }//loop over marked elements
524  }
525 
526 
527 protected:
528  virtual std::ostream& info(std::ostream& os) const{
529  return os << " AprioriVertex[ type = hp - Refinement Strategy ]";
530  //further output : history
531  // the parameters etc
532  }
533 
534 
535 private:
536 
537  //set of all vertice attributes at which h-refinement should be applied
539 
540  //current cells of given space
542 
543 };
544 
545 
546 }
547 
548 
549 
550 
551 
552 
553 
554 
555 
556 
557 
558 
559 
560 
561 
562 #endif // refinement_hh
virtual Connector & connector() const =0
Returns the connector.
void buildRefinement(const Marking &mark, const concepts::LocalEstimator< F > &estimator)
Definition: refinement.hh:494
virtual Connector0 * vertex(uint i) const =0
Returns a 0D component: vertex A 2D element consists of 1D and 0D elements, the vertices can be queri...
virtual Cell * child(uint i)=0
Returns a pointer to the ith child.
Abstract class for all refinement strategies in arbitray dimension.
Definition: refinement.hh:30
virtual std::ostream & info(std::ostream &os) const =0
Returns information in an output stream.
const Attribute & attrib() const
Returns the attribute of the connector.
Definition: connector.hh:108
concepts::Sequence< const concepts::Cell * > hRefCurCells_
Definition: refinement.hh:311
uint attrib() const
Returns the attribute.
Definition: connector.hh:38
AprioriVertex(const concepts::SpaceOnCells< F > &space, const uint attrb)
Definition: refinement.hh:473
#define conceptsException(exc)
Prepares an exception for throwing.
Definition: exceptions.hh:344
Prediction(const concepts::SpaceOnCells< F > &space, enum Mode mode=P, const Real hdec=sqrt(4.0), const Real pdec=sqrt(0.4), const Real nRdec=1.0)
Constructor for the prediction refinement strategy.
Definition: refinement.hh:122
virtual const ushort * innerDof(const concepts::Connector2 &cntr) const
Returns a pointer to the maximal polynomial degrees of the inner degree of freedom in the cell quad.
A cell in a mesh consist of topological information (neighbours, connectivity, orientation) and geome...
Definition: cell.hh:39
virtual std::ostream & info(std::ostream &os) const
Returns information in an output stream.
Definition: refinement.hh:282
virtual void apply(concepts::SpaceOnCells< F > &space)=0
Helper class for building 2D hp-FEM spaces (space pre builder).
concepts::HashMap< concepts::AdaptiveAdjustP< dim > > adjusts_
Definition: refinement.hh:87
void buildRefinement(const Marking &mark, const concepts::LocalEstimator< F > &estimator)
Definition: refinement.hh:187
#define conceptsAssert(cond, exc)
Assert that a certain condition is fulfilled.
Definition: exceptions.hh:394
virtual std::ostream & info(std::ostream &os) const
Returns information in an output stream.
Definition: refinement.hh:528
hp-adaptive space with 2D elements.
#define MAX_D
Definition: refinement.hh:19
concepts::HashMap< const concepts::Cell * > cells_
Definition: refinement.hh:306
concepts::Set< uint > attrbs_
Definition: refinement.hh:538
An array of objects.
Definition: bilinearForm.hh:23
Exception class for assertions.
Definition: exceptions.hh:258
const Key & key() const
Returns the key of the connector.
Definition: connector.hh:105
void rebuild(bool sameIndices=false)
Rebuilds the mesh and the elements due to adjustment orders.
A 2D element of the topology.
Definition: connector.hh:226
virtual bool vtxDof(const concepts::Connector2 &cntr, uint i) const
Returns, if the i th vertex of cell cntr has a degree of freedom.
Sequence with operations, output operator, and method of the particular element types.
Definition: sequence.hh:39
Abstract class on Marking.
Definition: marking.hh:15
virtual Scan * scan() const
Returns a scanner to iterate over the elements of the space.
concepts::HashMap< concepts::Array< ushort > > p_
Definition: refinement.hh:303
Abstract Class for h refinement only, i.e.
Definition: refinement.hh:59
const concepts::Set< uint > & marks() const
Definition: marking.hh:17
virtual const Cell & cell() const =0
Returns the cell on which the element is built.
Abstract class for a space.
Definition: space.hh:81
Class to describe adjustments to elements in an adaptive space.
Definition: hpMethod.hh:59
unsigned short ushort
Abbreviation for unsigned short.
Definition: typedefs.hh:48
Exception class to express a missing feature.
Definition: exceptions.hh:206
void addVertex(const uint attrb)
Definition: refinement.hh:489
Prediction hp - adaptive refinement strategy, proposed by Melenk and Wohlmuth.
Definition: refinement.hh:99
virtual Scanner * scan() const =0
Returns a scanner to iterate over the elements of the space.
virtual void apply(concepts::SpaceOnCells< F > &space)
Definition: refinement.hh:64
concepts::HashMap< const concepts::Cell * > cells_
Definition: refinement.hh:541
concepts::HashMap< Real > pred_
Definition: refinement.hh:300
uint key() const
Returns the key.
Definition: connector.hh:66
virtual void adjust(const concepts::Element< F > &elm, const concepts::AdaptiveAdjustP< 2 > &a)
Adjusts the space in the next rebuild step for this element.
virtual void apply(concepts::SpaceOnCells< F > &space)
Definition: refinement.hh:237
Class providing an output operator.
double Real
Type normally used for a floating point number.
Definition: typedefs.hh:36
A 0D element of the topology.
Definition: connector.hh:147
Page URL: http://wiki.math.ethz.ch/bin/view/Concepts/WebHome
21 August 2020
© 2020 Eidgenössische Technische Hochschule Zürich