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
52namespace 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
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 */
Definition femMisc.hpp:378
Definition femMisc.hpp:393
Definition femMisc.hpp:358
Definition femMisc.hpp:310
Definition femMisc.hpp:343
Definition femMisc.hpp:327
Definition femMisc.hpp:88
Definition femMisc.hpp:212
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
Definition femMisc.hpp:175
cvect(ccreal r)
cvect( cvect& r) { val[0]=r.val[0];val[1]=r.val[1];} // copy constructor
Definition femMisc.hpp:180

This is the FreeFEM reference manual
Provided by The KFEM project