FreeFem  3.5.x
femMisc.hpp
1 // Emacs will be in -*- Mode: c++ -*-
2 //
3 // ********** DO NOT REMOVE THIS BANNER **********
4 //
5 // SUMMARY: Language for a Finite Element Method
6 //
7 //
8 // AUTHORS: C. Prud'homme
9 // ORG :
10 // E-MAIL : prudhomm@users.sourceforge.net
11 //
12 // ORIG-DATE: June-94
13 // LAST-MOD: 12-Jul-01 at 09:50:55 by
14 //
15 // DESCRIPTION:
16 /*
17  This program is free software; you can redistribute it and/or modify
18  it under the terms of the GNU General Public License as published by
19  the Free Software Foundation; either version 2 of the License, or
20  (at your option) any later version.
21 
22  This program is distributed in the hope that it will be useful,
23  but WITHOUT ANY WARRANTY; without even the implied warranty of
24  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25  GNU General Public License for more details.
26 
27  You should have received a copy of the GNU General Public License
28  along with this program; if not, write to the Free Software
29  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
30 
31 */
32 // DESCRIP-END.
33 //
34 
35 // To change from comlex to real search everywhere for /* C2R */
36 
37 #ifndef __OPCLASS_H
38 #ifdef __GNUG__
39 #pragma interface
40 #endif
41 #define __OPCLASS_H
42 
43 #include <cmath>
44 #include <iostream>
45 
46 
47 #define mmax(a,b)(a>b?a:b)
48 #define mmin(a,b)(a<b?a:b)
49 #define ssqr(x) ((x)*(x))
50 
51 
52 namespace fem
53 {
54  extern int N;
55  extern int N2;
56 
57  class Complex;
58  class cvect;
59  class cmat;
60 
61 #define sizecmat sizeof(cmat)
62 #define sizecvect sizeof(cvect)
63 #define deter2(a11,a12,a21,a22) (a11*a22 - a12*a21)
64 
65  const float epsilon =(float)1.0e-20;
66 
67  //typedef float ccreal; typedef float creal;
68  typedef float ccreal;
69  typedef Complex creal; // Complex : change this (see also OPClass.c) /* C2R */
70 
71  typedef Complex ccomplex;
72 #define sizeccomplex sizeof(ccomplex)
73 
74 #define sizecreal sizeof(creal)
75 #define sizeccreal sizeof(ccreal)
76 
77  void cgauss( cmat& a, cvect& b);
78 
79  float norm2(const float& a);// { return a>0?a:-a; }
80  float imagpart(const float& a);//{ return 0; }
81  float realpart(const float& a);//{ return a; }
82  cmat id(const cvect& a);
83  Complex id(const Complex& x);
84  void myassert(int what);
85 
86  class Complex
87 
88  {
89  protected:
90  float re,im; // Internal variables
91  public:
92  Complex() { re = 0.F; im =0.F; } // default constructor
93  Complex(float r, float i =0.F) {re =r; im =i;} // copy constructor
94 
95  float& real () { return re;};
96  float real() const { return re;}
97  float& imag() { return im;};
98  float imag() const { return im;}
99 
100  Complex conjug() {return Complex(re,-im);}
101 
102  Complex& operator=(const Complex& a) { re =a.re; im = a.im; return *this; }
103  //Complex& operator=(const float& a) { re =a; im = 0.F; return *this; }
104  Complex& operator*=(const Complex& a) { re =a.re*re - a.im*im;im= a.re*im + a.im*re; return *this; }
105  Complex& operator/=(const Complex& a) {
106  float xl=ssqr(a.re) + ssqr(a.im); re =(a.re*re+a.im*im)/xl; im =(a.re*im-a.im*re)/xl;
107  return *this;
108  }
109  Complex& operator+=(const Complex& a) { re +=a.re; im += a.im; return *this; }
110  Complex& operator-=(const Complex& a) { re -=a.re; im -= a.im; return *this; }
111  Complex& operator*=(const float a) { re *= a;im *= a; return *this; }
112  Complex& operator/=(const float a) { re /= a;im /= a; return *this; }
113  Complex& operator+=(const float a) { re +=a; return *this; }
114  Complex& operator-=(const float a) { re -=a; return *this; }
115  Complex& operator-() { re =-re; im =-im; return *this; }
116  float modul2() const {return ssqr(re)+ssqr(im);}
117  float arg() const;
118 
119  friend Complex operator *(const Complex a, const Complex b) {return Complex(a.re*b.re - a.im*b.im, a.re*b.im + a.im*b.re);}
120  friend Complex operator *(const Complex a, const float b) { return Complex(a.re*b, a.im*b); }
121  friend Complex operator *(const float b, const Complex a) { return Complex(a.re*b, a.im*b); }
122  friend Complex operator /(const Complex a, const float b) { return Complex(a.re/b, a.im/b); }
123  friend Complex operator /(const float b, const Complex a) { return Complex(b,0.F)/a; }
124  friend Complex operator /(const Complex a, const Complex b) {
125  Complex c; float xl=ssqr(b.re) + ssqr(b.im); c.re =(a.re*b.re+a.im*b.im)/xl; c.im =(b.re*a.im-b.im*a.re)/xl;
126  return c;
127  }
128  friend Complex operator +(const Complex a, const Complex b) {return Complex(a.re+b.re, a.im+b.im);}
129  friend Complex operator +(const Complex a, const float b) { return Complex(a.re+b, a.im); }
130  friend Complex operator +(const float b, const Complex a) { return Complex(a.re+b, a.im); }
131  friend Complex operator -(const Complex a, const Complex b) {return Complex(a.re-b.re, a.im-b.im);}
132  friend Complex operator -(const Complex a, const float b) { return Complex(a.re-b, a.im); }
133  friend Complex operator -(const float b, const Complex a) { return Complex(b-a.re, -a.im); }
134  friend float norm2( const Complex& a) { return a.modul2(); }
135  friend float realpart(const Complex& a){ return a.re; }
136  friend float imagpart(const Complex& a){ return a.im; }
137  friend std::ostream& operator<<(std::ostream& os, const Complex& a);
138  friend std::istream& operator>>(std::istream& os, Complex& a);
139  friend Complex id(const Complex& x);
140 
141  friend Complex pow(const Complex&, const float&);
142  friend float cos(const Complex&);
143  friend float sin(const Complex&);
144 
145  };
146  inline float
147  Complex::arg() const
148  {
149  if (modul2() > 1e-8)
150  if (im > 0)
151  return acos(re/sqrt(modul2()));
152  else
153  return 8*atan(1.) - acos(re/sqrt(modul2()));
154  else
155  return 0.;
156 
157  }
158  inline Complex
159  pow(const Complex& c, const float& f)
160  {
161  Complex z(std::cos(f*c.arg()), std::sin(f*c.arg()));
162  return std::pow(std::sqrt(c.modul2()),f)*z;
163  }
164  inline float
165  cos(const Complex& c)
166  {
167  return std::cos(c.real())*std::cosh(c.imag()) + std::sin(c.real())*std::sinh(c.imag());\
168  }
169  inline float
170  sin(const Complex& c)
171  {
172  return std::sin(c.real())*std::cosh(c.imag()) - std::cos(c.real())*std::sinh(c.imag());
173  }
174  //____________________________________________________________________________
175  class cvect {public: ccreal val[2]; // Internal cvector to n values
176  //----------------------------------
177  cvect() { val[0]=0.F; val[1]=0.F; } // constructor
178  cvect(const cvect& r) { val[0]=r.val[0];val[1]=r.val[1];} // copy constructor
180  cvect(ccreal r) { val[0]=r;val[1]=r;} // copy constructor from a creal
181  ccreal& operator[] ( int i) { /*myassert((i< 2)&&(i>=0)); */ return val[i];}
182  const ccreal& operator[] ( int i) const { /*myassert((i< 2)&&(i>=0));*/ return val[i];}
183  float modul2() { return norm2(val[0])+norm2(val[1]);} // public fnct
184  ccreal operator*=(const cvect& a) { return val[0]*a.val[0] + val[1]*a.val[1]; }
185  cvect& operator*=(const ccreal a) { val[0] *= a; val[1] *= a; return *this; }
186  friend ccreal operator *(const cvect& a, const cvect& b) { cvect c(a); return c *= b; }
187  friend cvect operator *(const cvect& a, const ccreal b) { cvect c(a); return c *= b; }
188  friend cvect operator *(const ccreal b, const cvect& a) { cvect c(a); return c *= b; }
189  cvect& operator/=(const ccreal a) {val[0] /= a;val[1] /= a; return *this; }
190  friend cvect operator /(const cvect& a, const ccreal b) { cvect c(a); return c /= b; }
191  cvect& operator+=(const cvect& a) { val[0] += a.val[0];val[1] += a.val[1]; return *this; }
192  cvect& operator+=(const ccreal a) { val[0] += a;val[1] += a; return *this; }
193  friend cvect operator +(const cvect& a, const cvect& b) { cvect c(a); return c += b; }
194  friend cvect operator +(const cvect& a, const ccreal b) { cvect c(a); return c += b; }
195  friend cvect operator +(const ccreal b, const cvect& a) { cvect c(a); return c += b; }
196  cvect& operator-=(const cvect& a) { val[0] -= a.val[0];val[1] -= a.val[1]; return *this; }
197  cvect& operator-=(const ccreal a) { val[0] -= a;val[1] -= a; return *this; }
198  friend cvect operator -(const cvect& a, const cvect& b) { cvect c(a); return c -= b; }
199  friend cvect operator -(const cvect& a, const ccreal b) { cvect c(a); return c -= b; }
200  friend cvect operator -(const ccreal b, const cvect& a) { cvect c(b); return c -= a; }
201  cvect& operator=(const cvect& a) { val[0] = a.val[0];val[1] = a.val[1]; return *this; }
202  cvect& operator=(const ccreal a) {val[0] = a; val[1] = a; return *this; }
203  cvect& operator-() { val[0] = -val[0]; val[1] = -val[1]; return *this; }
204  friend float norm2( cvect& a) { return a.modul2();}
205  friend float realpart(const cvect& a){ return realpart(a.val[0]); }
206 
207  friend std::ostream& operator<<(std::ostream& os, cvect& a);
208  friend std::istream& operator>>(std::istream& os, cvect& a);
209  };
210 
211  //_______________________________
212  class cmat {
213  //---------
214  public: ccreal val[4]; // Internal cvector to n values
215  cmat() {val[0]=0.F;val[1]=0.F;val[2]=0.F;val[3]=0.F; } // constructor
216  cmat(const cmat& r) {val[0]=r.val[0];val[1]=r.val[1];val[2]=r.val[2];val[3]=r.val[3];} // copy constructor
218  cmat(ccreal r) { val[0] = r;val[1] = r;val[2] = r;val[3] = r;} // copy constructor from a creal produces a diag mat
219  ccreal& operator() (int i, int j) { /*myassert((i< N)&&(i>=0)&&(j<N)&&(j>=0));*/ return val[i+i+j];}
220  ccreal& operator[] (int i) { /*myassert((i< N2)&&(i>=0));*/ return val[i];}
221  const ccreal& operator[] (int i) const { /*myassert((i< N2)&&(i>=0));*/ return val[i];}
222  float modul2() {return val[0]+val[1]+val[2]+val[3];} // public fnct
223  cmat operator*=(const cmat& a)
224  { cmat b;
225  b.val[0] =val[0]*a.val[0]+val[1]*a.val[2];
226  b.val[1] =val[0]*a.val[1]+val[1]*a.val[3];
227  b.val[2] =val[2]*a.val[0]+val[3]*a.val[2];
228  b.val[3] =val[2]*a.val[1]+val[3]*a.val[3];
229  return b;
230  }
231  cmat& operator*=(const ccreal a) { val[0] *= a; val[1] *= a; val[2] *= a; val[3] *= a; return *this; }
232  friend cmat operator *(const cmat& a, const cmat& b) { cmat c(a); return c *= b; }
233  friend cmat operator *(const cmat& a, const ccreal b) { cmat c(a); return c *= b; }
234  friend cmat operator *(const ccreal b, const cmat& a) { cmat c(a); return c *= b; }
235  cmat& operator/=(const ccreal a) { val[0] /= a; val[1] /= a; val[2] /= a; val[3] /= a; return *this; }
236  friend cmat operator /(const cmat& a, const ccreal b) { cmat c(a); return c /= b; }
237  cmat& operator+=(const cmat& a){val[0]+=a.val[0];val[1]+=a.val[1];val[2]+=a.val[2];val[3]+=a.val[3]; return *this; }
238  cmat& operator+=(const ccreal a) { val[0] += a; val[1] += a; val[2] += a; val[3] += a; return *this; }
239  friend cmat operator +(const cmat& a, const cmat& b) { cmat c(a); return c += b; }
240  friend cmat operator +(const cmat& a, const ccreal b) { cmat c(a); return c += b; }
241  friend cmat operator +(const ccreal b, const cmat& a) { cmat c(a); return c += b; }
242 
243  cmat& operator-=(const cmat& a){val[0]-=a.val[0];val[1]-=a.val[1];val[2]-=a.val[2];val[3]-=a.val[3]; return *this; }
244  cmat& operator-=(const ccreal a) { val[0] -= a; val[2] -= a; val[1] -= a; val[3] -= a; return *this; }
245  friend cmat operator -(const cmat& a, const cmat& b) { cmat c(a); return c -= b; }
246  friend cmat operator -(const cmat& a, const ccreal b) { cmat c(a); return c -= b; }
247  friend cmat operator -(const ccreal b, const cmat& a) { cmat c(b); return c -= a; }
248  cmat& operator=(const cmat& a){val[0]=a.val[0];val[1]=a.val[1];val[2]=a.val[2]; val[3] = a.val[3]; return *this; }
249  cmat& operator=(const ccreal a) { val[0] = a; val[1] = a; val[2] = a; val[3] = a; return *this; }
250  cmat& operator-() { val[0] = -val[0];val[1] = -val[1];val[2] = -val[2];val[3] = -val[3]; return *this; }
251  friend float norm2( cmat& a) { return a.modul2();}
252 
253  friend cvect operator *(const cmat& a, const cvect& b)
254  { cvect c; c.val[0] = a.val[0]*b.val[0]+a.val[1]*b.val[1];
255  c.val[1] = a.val[2]*b.val[0]+a.val[3]*b.val[1];
256  return c ; }
257  friend cvect operator *(const cvect& b, const cmat& a)
258  { cvect c; c.val[0] = a.val[0]*b.val[0]+a.val[2]*b.val[1];
259  c.val[1] = a.val[1]*b.val[0]+a.val[3]*b.val[1];
260  return c ; } //transposed
261  friend cvect operator /(const cvect& b,const cmat& a)
262  { cvect c; ccreal s = deter2(a.val[0],a.val[1], a.val[2],a.val[3]);
263  if(norm2(s)<epsilon) { s = epsilon;}
264  c.val[0] = deter2(b.val[0],a.val[1], b.val[1],a.val[3]) / s;
265  c.val[1] = deter2(a.val[0],b.val[0], a.val[2],b.val[1]) / s;
266  return c;
267  }
268  friend cmat operator /(const cmat& b,const cmat& a)
269  {
270  cmat c; ccreal s = deter2(a.val[0],a.val[1], a.val[2],a.val[3]);
271  if(norm2(s)<epsilon){ s = epsilon;}
272  c.val[0] = deter2(b.val[0],a.val[1], b.val[2],a.val[3]) / s;
273  c.val[2] = deter2(a.val[0],b.val[0], a.val[2],b.val[2]) / s;
274  c.val[1] = deter2(b.val[1],a.val[1], b.val[3],a.val[3]) / s;
275  c.val[3] = deter2(a.val[0],b.val[1], a.val[2],b.val[3]) / s;
276  return c;
277  }
278  friend std::ostream& operator<<(std::ostream& os, cmat& a);
279  friend std::istream& operator>>(std::istream& os, cmat& a);
280  friend cmat id(const cvect& a);
281  };
282 
283  /*
284  class Afloat
285  { public:
286  int szz;
287  float* cc;
288 
289  Afloat():szz(0),cc(NULL){}
290  Afloat(int sz=0) {cc = 0;
291  if(sz>0){ cc = new float[sz];
292  if(!cc) erreur("Out of Memory");
293  for(int i=0;i<sz;i++)cc[i]=0; }
294  szz = sz;
295  }
296  Afloat(Afloat& a) { cc = 0; szz = 0;
297  if(a.szz>0) { szz = a.szz; cc = new float[szz];
298  if(!cc) erreur("Out of Memory");
299  else for(int i=0;i<szz;i++)cc[i]=a.cc[i];}
300  }
301  ~Afloat() { delete [] cc;cc=0;szz = 0;}
302  float& operator[] (int i) { myassert((i< szz)&&(i>=0)); return cc[i];}
303  void init(int newSize) { myassert(!(szz || cc)); szz =newSize;
304  cc =new float[szz]; if(!cc) erreur("Out of Memory");
305  for(int i=0;i<szz;i++)cc[i]=0;;
306  }
307  };
308  */
309  class Acreal
310  {
311  public:
312  long szz;
313  creal* cc;
314 
315  Acreal(long=0);
316  Acreal(const Acreal&);
317  ~Acreal() { delete [] cc;cc=0;szz = 0;}
318  void destroy() {delete [] cc;cc=0;szz = 0;}
319  creal& operator[] (long i) { /*myassert((i< szz)&&(i>=0));*/ return cc[i];}
320  creal* operator&() { return cc;}
321  void init(long newSize);/* { myassert(!(szz || cc)); szz =newSize;
322  cc =new creal[szz]; if(!cc) erreur("Out of Memory");
323  for(int i=0;i<szz;i++)cc[i]=0;
324  } */
325  };
326  class Aint
327  {
328  public:
329  long szz;
330  int* cc;
331 
332  Aint(long=0);
333  Aint(const Aint&);
334  ~Aint() { delete [] cc;cc=0;szz = 0;}
335  void destroy() {delete [] cc;cc=0;szz = 0;}
336  int& operator[] (long i) { myassert((i< szz)&&(i>=0)); return cc[i];}
337  int* operator&() { return cc;}
338  void init (long);
339  };
340 
341 
342  class Acvect
343  {
344  public:
345  long szz;
346  cvect* cc;
347 
348  Acvect(long=0);
349  Acvect(const Acvect&);
350  ~Acvect() { delete [] cc;cc=0;szz = 0;}
351  void destroy() { delete [] cc;cc=0;szz = 0;}
352  cvect& operator[] (long i) { /*myassert((i< szz)&&(i>=0));*/ return cc[i];}
353  cvect* operator&() { return cc;}
354  void init (long);
355  };
356 
357  class Acmat
358  {
359  public:
360  long szz;
361  cmat* cc;
362 
363  Acmat(long=0);
364  Acmat(const Acmat&);
365  ~Acmat()
366  {
367  delete [] cc;
368  cc=0;
369  szz = 0;
370  }
371  void destroy() {delete [] cc;cc=0;szz = 0;}
372  cmat& operator[] (long i) { /*myassert((i< szz)&&(i>=0));*/ return cc[i];}
373  cmat* operator&() { return cc;}
374  void init (long);
375  };
376 
377  class AAcmat
378  {
379  public:
380  long szz;
381  Acmat* cc;
382 
383  AAcmat(long=0);
384  AAcmat(const AAcmat&);
385  ~AAcmat() { delete [] cc;cc=0;szz = 0;}
386  void destroy() {delete [] cc;cc=0;szz = 0;}
387  Acmat& operator[] (long i) { /*myassert((i< szz)&&(i>=0));*/ return cc[i];}
388  Acmat* operator&() { return cc;}
389  void init(long );
390  };
391 
392  class AAcreal
393  {
394  public:
395  long szz;
396  Acreal* cc;
397 
398  AAcreal(long=0);
399  AAcreal(const AAcreal& );
400  ~AAcreal() { delete [] cc;cc=0;szz = 0;}
401  void destroy() {delete [] cc;cc=0;szz = 0;}
402  Acreal& operator[] (long i) { /*myassert((i< szz)&&(i>=0));*/ return cc[i];}
403  Acreal* operator&() { return cc;}
404  void init (long);
405  };
406 }
407 
408 #endif /* __OPCLASS_H */
fem::Acvect
Definition: femMisc.hpp:342
fem::AAcreal
Definition: femMisc.hpp:392
fem::AAcmat
Definition: femMisc.hpp:377
fem::Complex
Definition: femMisc.hpp:86
fem::Aint
Definition: femMisc.hpp:326
fem::cmat::cmat
cmat(ccreal r)
cmat(cmat& r) { val[0]=r.val[0]; val[1]=r.val[1]; val[2]=r.val[2]; val[3]=r.val[3]; } // copy constru...
Definition: femMisc.hpp:218
fem::Acreal
Definition: femMisc.hpp:309
fem::cvect
Definition: femMisc.hpp:175
fem::cmat
Definition: femMisc.hpp:212
fem::cvect::cvect
cvect(ccreal r)
cvect( cvect& r) { val[0]=r.val[0];val[1]=r.val[1];} // copy constructor
Definition: femMisc.hpp:180
fem::Acmat
Definition: femMisc.hpp:357

This is the FreeFEM reference manual
Provided by The KFEM project