1 // ********** DO NOT REMOVE THIS BANNER **********
2 // ORIG-DATE: Jan 2008
3 // -*- Mode : c++ -*-
4 //
5 // SUMMARY : Generic Quadrature Formular
6 // USAGE : LGPL
7 // ORG : LJLL Universite Pierre et Marie Curi, Paris, FRANCE
8 // AUTHOR : Frederic Hecht
9 // E-MAIL : frederic.hecht@ann.jussieu.fr
10 //
11
12 /*
13
14 This file is part of Freefem++
15
16 Freefem++ is free software; you can redistribute it and/or modify
17 it under the terms of the GNU Lesser General Public License as published by
18 the Free Software Foundation; either version 2.1 of the License, or
19 (at your option) any later version.
20
21 Freefem++ is distributed in the hope that it will be useful,
22 but WITHOUT ANY WARRANTY; without even the implied warranty of
23 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 GNU Lesser General Public License for more details.
25
26 You should have received a copy of the GNU Lesser General Public License
27 along with Freefem++; if not, write to the Free Software
28 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
29
30 Thank to the ARN FF2A3 grant
31 ref:ANR-07-CIS7-002-01
32 */
33
34
35 #ifndef _QuadratureFormular_h
36 #define _QuadratureFormular_h
37
38 #include <iostream>
39
40 namespace Fem2D {
41 using namespace std;
42 #include "R3.hpp"
43
44
45
46 //class QuadratureFormular;
47 struct QuadratureWeight {
48 R a;
QuadratureWeightFem2D::QuadratureWeight49 QuadratureWeight(R aa): a(aa){}
50 };
51
52 template<class Rd>
53 class GQuadraturePoint: public QuadratureWeight,public Rd {
54 public:
55 typedef GQuadraturePoint QP ;
GQuadraturePoint()56 GQuadraturePoint(): QuadratureWeight(0),Rd() {}
GQuadraturePoint(R aa,const Rd & xx)57 GQuadraturePoint(R aa, const Rd &xx): QuadratureWeight(aa),Rd(xx) {}
GQuadraturePoint(const Rd & xx,R aa)58 GQuadraturePoint(const Rd & xx,R aa): QuadratureWeight(aa),Rd(xx) {}
operator R() const59 operator R() const {return a;}
GQuadraturePoint(R aa,R xx)60 GQuadraturePoint(R aa,R xx):QuadratureWeight(aa),Rd(xx) {}
GQuadraturePoint(R aa,R x,R y)61 GQuadraturePoint(R aa,R x,R y):QuadratureWeight(aa),Rd(x,y) {}
GQuadraturePoint(R aa,R x,R y,R z)62 GQuadraturePoint(R aa,R x,R y,R z):QuadratureWeight(aa),Rd(x,y,z) {}
63 template<class RD >
64 GQuadraturePoint<RD>
Bary(RD * K,double mes)65 Bary(RD * K, double mes) { return GQuadraturePoint<RD>(this->Rd::Bary(K),a*mes);}
66 };
67
68 template<class Rdd>
69 class GQuadratureFormular {
70
71 public:
72 typedef Rdd Rd;
73 typedef GQuadraturePoint<Rd> QuadraturePoint;
74 typedef GQuadraturePoint<Rd> QP;
75 int exact; // exact
76 int n; // nombre de point d'integration
77 const int size; // size of the array
78 private:
79 QP *p; // les point d'integration
80 const bool clean;
81 public:
82 // -- les fonctions ------------------
83 void Verification(); // for verification
GQuadratureFormular(int e,int NbOfNodes,QuadraturePoint * pp,bool c=false)84 GQuadratureFormular (int e,int NbOfNodes,QuadraturePoint *pp,bool c=false)
85 :exact(e), n(NbOfNodes),size(n),p(pp),clean(c) {Verification();}
GQuadratureFormular(int e,int NbOfNodes,const QuadraturePoint * pp,bool c=false)86 GQuadratureFormular (int e,int NbOfNodes,const QuadraturePoint *pp,bool c=false)
87 :exact(e),n(NbOfNodes),p(pp),clean(c) {Verification();}
88
GQuadratureFormular(int ex,QP p0,QP p1,QP p2,QP p3,QP p4)89 GQuadratureFormular(int ex,QP p0,QP p1,QP p2,QP p3,QP p4)
90 : exact(ex),n(5),size(n),p(new QP[5]),clean(true) { p[0]=p0;p[1]=p1;p[2]=p2;p[3]=p3;p[4]=p4;Verification();}
GQuadratureFormular(int ex,QP p0,QP p1,QP p2,QP p3)91 GQuadratureFormular(int ex,QP p0,QP p1,QP p2,QP p3)
92 : exact(ex),n(4),size(n),p(new QP[4]) ,clean(true){ p[0]=p0,p[1]=p1,p[2]=p2;p[3]=p3;Verification();}
GQuadratureFormular(int ex,QP p0,QP p1,QP p2)93 GQuadratureFormular(int ex,QP p0,QP p1,QP p2)
94 : exact(ex),n(3),size(n),p(new QP[3]),clean(true) { p[0]=p0,p[1]=p1,p[2]=p2;Verification();}
GQuadratureFormular(int ex,QP p0,QP p1)95 GQuadratureFormular(int ex,QP p0,QP p1)
96 : exact(ex),n(2),size(n),p(new QP[2]),clean(true) { p[0]=p0,p[1]=p1;Verification();}
GQuadratureFormular(int ex,QP p0)97 GQuadratureFormular(int ex,QP p0)
98 : exact(ex),n(1),size(n),p(new QP[1]),clean(true) { p[0]=p0;Verification();}
99 // bluid a empty GQuadratureFormular
GQuadratureFormular(int ssize)100 GQuadratureFormular(int ssize):exact(0),n(0),size(ssize),p(new QP[size]),clean(true) {}
101
operator [](int i) const102 const QP & operator [](int i) const {return p[i];}
operator ()(int i) const103 const QP & operator ()(int i) const {return p[i];}
~GQuadratureFormular()104 ~GQuadratureFormular() {if(clean) delete [] p;}
105
GQuadratureFormular(const GQuadratureFormular & QF,int mul=1)106 GQuadratureFormular(const GQuadratureFormular & QF, int mul=1)
107 :exact(QF.exact),n(QF.n),size(QF.size*mul),p(new QP[size*mul]),clean(true){ operator=(QF);}
operator =(const GQuadratureFormular & QF)108 void operator=(const GQuadratureFormular &QF)
109 {
110 assert(size>=QF.n);
111 n = QF.n;
112 for(int i=0;i<n;++i)
113 p[i]=QF.p[i];
114 }
operator *=(double c)115 void operator*=( double c)
116 {
117 for(int i=0;i<n;++i)
118 p[i].a *=c;
119 }
120 // Add new GQuadratureFormular on element K to the current Quadarture formular ..
121 // FH april 2014 ..
122 // to compute int under levelset ..
reset()123 void reset() { n =0;exact=0;}
124 template<class RD >
AddQFK(const GQuadratureFormular<RD> & QF,Rd * K,double mes,int n0=0)125 void AddQFK(const GQuadratureFormular<RD> &QF ,Rd *K,double mes,int n0=0)
126 {
127
128 ffassert( size >= n0 + QF.n );
129 n0 += n;
130 n = n0 + QF.n;
131 for(int i=0;i<QF.n;++i)
132 p[i+n0]=QF.p[i].Bary(K,mes);
133
134 }
135 private:
136 /* GQuadratureFormular(const GQuadratureFormular &)
137 :exact(0),n(0),p(0){assert(0);}
138 void operator=(const GQuadratureFormular &)
139 {assert(0);}
140 GQuadratureFormular()
141 :exact(0),n(0),p(0){assert(0);}*/
142 static const GQuadratureFormular * Default;
143 };
144
145
146
147 template<class Rd>
148 ostream& operator <<(ostream& , const GQuadratureFormular<Rd> & ) ;
149 template<class Rd>
150 ostream& operator <<(ostream& , GQuadraturePoint<Rd> & );
151
152 typedef GQuadratureFormular<R1> QuadratureFormular1d;
153
154
155 extern const QuadratureFormular1d QF_GaussLegendre1;
156 extern const QuadratureFormular1d QF_GaussLegendre2;
157 extern const QuadratureFormular1d QF_GaussLegendre3;
158 extern const QuadratureFormular1d QF_GaussLegendre4;
159 extern const QuadratureFormular1d QF_GaussLegendre5;
160 extern const QuadratureFormular1d QF_LumpP1_1D;
161
162
163 extern const GQuadratureFormular<R2> QuadratureFormular_T_1;
164 extern const GQuadratureFormular<R2> QuadratureFormular_T_1lump;
165 extern const GQuadratureFormular<R2> QuadratureFormular_T_2;
166 extern const GQuadratureFormular<R2> QuadratureFormular_T_5;
167 extern const GQuadratureFormular<R2> QuadratureFormular_T_2_4P1;
168 extern const GQuadratureFormular<R2> QuadratureFormular_T_7;
169 extern const GQuadratureFormular<R2> QuadratureFormular_T_9;
170
171 extern const GQuadratureFormular<R3> QuadratureFormular_Tet_1;
172 extern const GQuadratureFormular<R3> QuadratureFormular_Tet_1lump;
173 extern const GQuadratureFormular<R3> QuadratureFormular_Tet_2;
174 extern const GQuadratureFormular<R3> QuadratureFormular_Tet_5;
175
176
177
178 template<class Rd>
179 GQuadratureFormular<Rd> * QF_Simplex(int exact);
180 // { return QF_exact<GQuadratureFormular<Rd>,Rd::d+1>(exact);}
181
182 template<class Rd>
setQF(GQuadratureFormular<Rd> & FI,const GQuadratureFormular<Rd> & FI1,const GQuadratureFormular<Rd> & FI0,Rd Q[Rd::d][Rd::d+1],double * cmes,int n)183 void setQF( GQuadratureFormular<Rd> &FI,
184 const GQuadratureFormular<Rd> & FI1 ,
185 const GQuadratureFormular<Rd> & FI0 ,
186 Rd Q[Rd::d][Rd::d+1],
187 double *cmes,
188 int n)
189 {
190 FI.reset();
191 for(int i=0; i< n; ++i)
192 if(cmes[i]==1)
193 FI=FI1;
194 else if(cmes[i] > 1e-4)
195 FI.AddQFK(FI1,Q[i],cmes[i]);
196 else if ( cmes[i]> 1e-8 )
197 FI.AddQFK(FI0,Q[i],cmes[i]);
198
199 }
200
201
202
203 }
204
205 namespace Fem2D {
206 typedef GQuadratureFormular<R2> QuadratureFormular;
207 typedef GQuadraturePoint<R2> QuadraturePoint;
208 typedef GQuadratureFormular<R1> QuadratureFormular1d;
209 typedef GQuadraturePoint<R1> QuadratureFormular1dPoint;
210 GQuadraturePoint<R1> * GaussLegendre(int nn);
211
212 }
213 #endif
214