LORENE
mg3d.C
1 /*
2  * Methods of class Mg3d
3  *
4  */
5 
6 /*
7  * Copyright (c) 1999-2000 Jean-Alain Marck
8  * Copyright (c) 1999-2001 Eric Gourgoulhon
9  *
10  * This file is part of LORENE.
11  *
12  * LORENE is free software; you can redistribute it and/or modify
13  * it under the terms of the GNU General Public License as published by
14  * the Free Software Foundation; either version 2 of the License, or
15  * (at your option) any later version.
16  *
17  * LORENE is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  * GNU General Public License for more details.
21  *
22  * You should have received a copy of the GNU General Public License
23  * along with LORENE; if not, write to the Free Software
24  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25  *
26  */
27 
28 
29 
30 
31 /*
32  * $Id: mg3d.C,v 1.22 2023/05/24 09:48:50 g_servignat Exp $
33  * $Log: mg3d.C,v $
34  * Revision 1.22 2023/05/24 09:48:50 g_servignat
35  * Added plus_half grid in angular direction for dealiasing
36  *
37  * Revision 1.21 2020/01/27 11:00:19 j_novak
38  * New include <stdexcept> to be compatible with older versions of g++
39  *
40  * Revision 1.20 2018/12/05 15:03:20 j_novak
41  * New Mg3d constructor from a formatted file.
42  *
43  * Revision 1.19 2016/12/05 16:17:59 j_novak
44  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
45  *
46  * Revision 1.18 2014/10/13 08:53:07 j_novak
47  * Lorene classes and functions now belong to the namespace Lorene.
48  *
49  * Revision 1.17 2014/10/06 15:13:14 j_novak
50  * Modified #include directives to use c++ syntax.
51  *
52  * Revision 1.16 2013/06/05 15:00:26 j_novak
53  * Suppression of all classes derived from Grille3d. Now Grille3d is no
54  * longer an abstract class. r-samplings are only one of RARE, FIN or
55  * UNSURR (FINJAC has been removed). Instead, Mg3d possesses a new member
56  * colloc_r[nzone] defining the type of collocation points (spectral
57  * bases) in each domain.
58  *
59  * Revision 1.15 2012/01/17 10:37:42 j_penner
60  * added a constructor that treats all domains as type FIN
61  *
62  * Revision 1.14 2008/08/19 06:42:00 j_novak
63  * Minor modifications to avoid warnings with gcc 4.3. Most of them concern
64  * cast-type operations, and constant strings that must be defined as const char*
65  *
66  * Revision 1.13 2007/12/11 15:28:15 jl_cornou
67  * Jacobi(0,2) polynomials partially implemented
68  *
69  * Revision 1.12 2006/05/17 13:17:03 j_novak
70  * New member g_angu_1dom, the one-domain angular grid associated with the
71  * current grid.
72  *
73  * Revision 1.11 2005/10/07 08:47:21 j_novak
74  * Addition of the pointer g_non_axi on a grid, with at least 5 points in the
75  * theta direction and 4 in the phi one (for tensor rotations).
76  *
77  * Revision 1.10 2004/07/06 13:36:29 j_novak
78  * Added methods for desaliased product (operator |) only in r direction.
79  *
80  * Revision 1.9 2003/10/27 16:21:54 e_gourgoulhon
81  * Treated the special case nz=1 in the simplified constructor.
82  *
83  * Revision 1.8 2003/06/20 14:50:15 f_limousin
84  * Add the operator==
85  *
86  * Revision 1.7 2003/06/18 08:45:27 j_novak
87  * In class Mg3d: added the member get_radial, returning only a radial grid
88  * For dAlembert solver: the way the coefficients of the operator are defined has been changed.
89  *
90  * Revision 1.6 2002/10/16 14:36:42 j_novak
91  * Reorganization of #include instructions of standard C++, in order to
92  * use experimental version 3 of gcc.
93  *
94  * Revision 1.5 2002/05/07 07:36:03 e_gourgoulhon
95  * Compatibilty with xlC compiler on IBM SP2:
96  * suppressed the parentheses around argument of instruction new:
97  * e.g. t = new (Tbl *[nzone]) --> t = new Tbl*[nzone]
98  *
99  * Revision 1.4 2001/12/12 09:23:46 e_gourgoulhon
100  * Parameter compact added to the simplified constructor of class Mg3d
101  *
102  * Revision 1.3 2001/12/11 06:48:30 e_gourgoulhon
103  * Addition of the simplified constructor
104  *
105  * Revision 1.2 2001/12/04 21:27:54 e_gourgoulhon
106  *
107  * All writing/reading to a binary file are now performed according to
108  * the big endian convention, whatever the system is big endian or
109  * small endian, thanks to the functions fwrite_be and fread_be
110  *
111  * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
112  * LORENE
113  *
114  * Revision 2.10 2001/05/26 14:50:46 eric
115  * *** empty log message ***
116  *
117  * Revision 2.9 2001/05/26 13:25:59 eric
118  * Ajout du membre g_twice (grille double pour le desaliasing)
119  * Modif de la declaration de g_angu (pointeur mutable)
120  * g_twice et g_angu ne sont calcules que si necessaire (cad si
121  * on appelle la fonction get_twice() ou get_angu()).
122  *
123  * Revision 2.8 2000/03/22 13:38:51 eric
124  * Remplacement des iendl par endl dans <<
125  *
126  * Revision 2.7 1999/10/12 15:04:29 eric
127  * *** empty log message ***
128  *
129  * Revision 2.6 1999/10/12 15:03:30 eric
130  * *** empty log message ***
131  *
132  * Revision 2.5 1999/09/30 14:58:16 eric
133  * Operator!= declare const
134  *
135  * Revision 2.4 1999/09/30 14:12:04 eric
136  * sauve declaree const.
137  *
138  * Revision 2.3 1999/09/30 12:52:52 eric
139  * Depoussierage.
140  * Documentation.
141  *
142  * Revision 2.2 1999/03/01 14:35:21 eric
143  * Modif affichage (operator<<)
144  *
145  *
146  * $Header: /cvsroot/Lorene/C++/Source/Mg3d/mg3d.C,v 1.22 2023/05/24 09:48:50 g_servignat Exp $
147  *
148  */
149 // C++ Headers
150 #include <stdexcept>
151 
152 // C Headers
153 #include <cstdlib>
154 #include <cassert>
155 
156 // Lorene headers
157 #include "grilles.h"
158 #include "type_parite.h"
159 #include "utilitaires.h"
160 
161  //--------------//
162  // Multi-grille //
163  //--------------//
164 
165 
166 //=============================================================================
167 // General constructor
168 //=============================================================================
169 
170 
171 namespace Lorene {
172 Mg3d::Mg3d(int nz, int nbr[], int typr[], int nbt[], int typt, int nbp[],
173  int typp, int* base_r)
174  : nzone(nz), type_t(typt), type_p(typp)
175 {
176 
177  // Type d'echantillonnage dans chaque zone
178  type_r = new int[nz];
179  colloc_r = new int[nz] ;
180  bool cheb = (base_r == 0x0) ;
181  for (int i=0 ; i<nz ; i++) {
182  type_r[i] = typr[i];
183  colloc_r[i] = cheb ? BASE_CHEB : base_r[i] ;
184  }
185 
186  // Nombre de points
187  nr = new int[nz];
188  nt = new int[nz];
189  np = new int[nz];
190  for (int i=0 ; i<nz ; i++) {
191  nr[i] = nbr[i] ;
192  nt[i] = nbt[i] ;
193  np[i] = nbp[i] ;
194  }
195 
196  // Les grilles
197  // -----------
198  g = new Grille3d*[nz] ;
199 
200  for (int i=0; i<nz; i++) {
201 
202  g[i] = new Grille3d(nr[i], nt[i], np[i], type_r[i], type_t, type_p,
203  colloc_r[i]) ;
204  } // fin de la boucle sur les zones
205 
206  // Pointers on derived grids initiated to 0x0:
207  // -------------------------------------------
208 
209  set_deriv_0x0() ;
210 
211 
212 }
213 
214 //=============================================================================
215 // Simplified constructor
216 //=============================================================================
217 
218 Mg3d::Mg3d(int nz, int nbr, int nbt, int nbp, int typt, int typp,
219  bool compact, bool leg)
220  : nzone(nz),
221  type_t(typt),
222  type_p(typp) {
223 
224  // Type of r sampling in each domain:
225  type_r = new int[nz];
226  colloc_r = new int[nz];
227  type_r[0] = RARE ;
228  colloc_r[0] = leg ? BASE_LEG : BASE_CHEB ;
229  for (int l=1; l<nz-1; l++) {
230  type_r[l] = FIN ;
231  colloc_r[l] = leg ? BASE_LEG : BASE_CHEB ;
232  }
233  if (nz > 1) {
234  if (compact) {
235  type_r[nz-1] = UNSURR ;
236  colloc_r[nz-1] = BASE_CHEB ;
237  }
238  else {
239  type_r[nz-1] = FIN ;
240  colloc_r[nz-1] = leg ? BASE_LEG : BASE_CHEB ;
241  }
242  }
243 
244  // Same number of points in all domains:
245  nr = new int[nz];
246  nt = new int[nz];
247  np = new int[nz];
248  for (int l=0 ; l<nz ; l++) {
249  nr[l] = nbr ;
250  nt[l] = nbt ;
251  np[l] = nbp ;
252  }
253 
254  // Les grilles
255  // -----------
256  g = new Grille3d*[nz] ;
257 
258  for (int i=0; i<nz; i++) {
259 
260  g[i] = new Grille3d(nr[i], nt[i], np[i], type_r[i], type_t, type_p,
261  colloc_r[i]) ;
262  } // fin de la boucle sur les zones
263 
264  // Pointers on derived grids initiated to 0x0:
265  // -------------------------------------------
266 
267  set_deriv_0x0() ;
268 
269 }
270 
271 //=============================================================================
272 // Simplified shell constructor
273 // Note: This does not handle the nucleus or the CED!
274 //=============================================================================
275 
276 Mg3d::Mg3d(int nz, int nbr, int nbt, int nbp, int typt, int typp)
277  : nzone(nz),
278  type_t(typt),
279  type_p(typp) {
280 
281  // Type of r sampling in each domain:
282  type_r = new int[nz];
283  colloc_r = new int[nz] ;
284  for (int l=0; l<nz; l++) {
285  type_r[l] = FIN ;
286  colloc_r[l] = BASE_CHEB ;
287  }
288 
289  // Same number of points in all domains:
290  nr = new int[nz];
291  nt = new int[nz];
292  np = new int[nz];
293  for (int l=0 ; l<nz ; l++) {
294  nr[l] = nbr ;
295  nt[l] = nbt ;
296  np[l] = nbp ;
297  }
298 
299  // Les grilles
300  // -----------
301  g = new Grille3d*[nz] ;
302 
303  for (int i=0; i<nz; i++) {
304 
305  g[i] = new Grille3d(nr[i], nt[i], np[i], type_r[i], type_t, type_p,
306  colloc_r[i]) ;
307 
308  } // fin de la boucle sur les zones
309 
310  // Pointers on derived grids initiated to 0x0:
311  // -------------------------------------------
312 
313  set_deriv_0x0() ;
314 
315 }
316 
317 //=============================================================================
318 // Constructors from a file
319 //=============================================================================
320 
321 
322  // From a formatted file...
323  //=========================
324  Mg3d::Mg3d(const string& filename) {
325 
326  // Opening of the file
327  //-----------------------------------------------------------------------
328  ifstream parfile(filename.c_str()) ;
329  if (!parfile) {
330  string message = "Unable to open file " ;
331  message += filename ;
332  throw ios_base::failure(message);
333  };
334 
335  string tmp_str = "Definition of the Mg3d" ;
336  if (!search_file(parfile, tmp_str)) {
337  string mesg = "No data found to contruct an object Mg3d in " ;
338  mesg += filename ;
339  throw invalid_argument(mesg) ;
340  }
341  // parfile >> tmp_str ;
342  //cout << tmp_str ;
343  parfile.ignore(1000, '\n') ;
344 
345  // Number of domains
346  //---------------------
347  parfile >> nzone ; parfile.ignore(1000, '\n') ; // Number of domains
348  parfile >> tmp_str ;
349 
350  // Theta & phi symmetries + number of grid points
351  //-----------------------------------------------
352  if ( tmp_str.compare("SYM") == 0) {
353  type_t = SYM ;
354  }
355  else if ( tmp_str.compare("NONSYM") == 0)
356  {
357  type_t = NONSYM ;
358  }
359  else
360  throw
361  invalid_argument("Mg3d::Mg3d(string): incorrect value for theta symmetry") ;
362  parfile.ignore(1000, '\n') ;
363 
364  parfile >> tmp_str ;
365  if ( tmp_str.compare("SYM") == 0)
366  {
367  type_p = SYM ;
368  }
369  else if ( tmp_str.compare("NONSYM") == 0)
370  {
371  type_p = NONSYM ;
372  }
373  else
374  throw
375  invalid_argument("Mg3d::Mg3d(string): incorrect value for phi symmetry") ;
376  parfile.ignore(1000, '\n') ;
377 
378  int nbp ; parfile >> nbp ; parfile.ignore(1000, '\n') ;
379  int nbt ; parfile >> nbt ; parfile.ignore(1000, '\n') ;
380 
381  // Radial number of points and sampling type
382  //---------------------------------------------
383  nr = new int[nzone] ;
384  nt = new int[nzone] ;
385  np = new int[nzone] ;
386  type_r = new int[nzone] ;
387 
388  for (int i=0; i<nzone; i++) {
389  parfile >> nr[i] >> tmp_str ;
390  assert (nr[i] > 0) ;
391  if ( tmp_str.compare("nucleus") == 0)
392  {
393  type_r[i] = RARE ;
394  }
395  else if ( tmp_str.compare("shell") == 0)
396  {
397  type_r[i] = FIN ;
398  }
399  else if ( tmp_str.compare("ced") == 0)
400  {
401  type_r[i] = UNSURR ;
402  }
403  else
404  throw
405  invalid_argument("Mg3d::Mg3d(string): incorrect value for the type of domain") ;
406  parfile.ignore(1000, '\n') ;
407  nt[i] = nbt ;
408  np[i] = nbp ;
409  } // end of loop on domains
410 
411  colloc_r = new int[nzone] ;
412  bool legendre ;
413  parfile >> legendre ;
414  for (int i=0; i<nzone; i++)
415  colloc_r[i] = (legendre ? BASE_LEG : BASE_CHEB) ;
416 
417 
418  // Grids
419  // -----------
420  g = new Grille3d*[nzone] ;
421 
422  for (int i=0; i<nzone; i++) {
423 
424  g[i] = new Grille3d(nr[i], nt[i], np[i], type_r[i], type_t, type_p,
425  colloc_r[i]) ;
426 
427  } // end of loop on domains
428 
429  // Pointers on derived grids initiated to 0x0:
430  // -------------------------------------------
431 
432  set_deriv_0x0() ;
433 
434  }
435 
436 
437 /*
438  * Construction a partir d'un fichier.
439  * Cette facon de faire est abominable. Cependant je ne vois pas comment
440  * faire autrement... j.a.
441  */
442 Mg3d::Mg3d(FILE* fd, bool read_base)
443 {
444  // Lecture sur le fichier
445  fread_be(&nzone, sizeof(int), 1, fd) ; // nzone
446  nr = new int[nzone] ;
447  fread_be(nr, sizeof(int), nzone, fd) ; // nr
448  nt = new int[nzone] ;
449  fread_be(nt, sizeof(int), nzone, fd) ; // nt
450  np = new int[nzone] ;
451  fread_be(np, sizeof(int), nzone, fd) ; // np
452  type_r = new int[nzone] ;
453  fread_be(type_r, sizeof(int), nzone, fd) ; // type_r
454  fread_be(&type_t, sizeof(int), 1, fd) ; // type_t
455  fread_be(&type_p, sizeof(int), 1, fd) ; // type_p
456  colloc_r = new int[nzone] ;
457  if (read_base)
458  fread_be(colloc_r, sizeof(int), nzone, fd) ; // colloc_r
459 
460  // Les grilles
461  // -----------
462 
463  g = new Grille3d*[nzone] ;
464  for (int i=0; i<nzone; i++) {
465  if (!read_base) colloc_r[i] = BASE_CHEB ;
466  g[i] = new Grille3d(nr[i], nt[i], np[i], type_r[i], type_t, type_p,
467  colloc_r[i]) ;
468 
469  } // fin de la boucle sur les zones
470 
471  // Pointers on derived grids initiated to 0x0:
472  // -------------------------------------------
473 
474  set_deriv_0x0() ;
475 
476 }
477 
478 // Destructeur
479 // -----------
481 
482  del_deriv() ; // Deletes the derived quantities
483 
484  delete [] nr ;
485  delete [] nt ;
486  delete [] np ;
487  delete [] type_r ;
488  delete [] colloc_r ;
489  for (int i=0 ; i<nzone ; i++) {
490  delete g[i] ;
491  }
492  delete [] g ;
493 
494 }
495 
496 //==================================================================
497 // Write in a file
498 //==================================================================
499 
500 void Mg3d::sauve(FILE* fd, bool save_base) const {
501  fwrite_be(&nzone, sizeof(int), 1, fd) ; // nzone
502  fwrite_be(nr, sizeof(int), nzone, fd) ; // nr
503  fwrite_be(nt, sizeof(int), nzone, fd) ; // nt
504  fwrite_be(np, sizeof(int), nzone, fd) ; // np
505  fwrite_be(type_r, sizeof(int), nzone, fd) ; // type_r
506  fwrite_be(&type_t, sizeof(int), 1, fd) ; // type_t
507  fwrite_be(&type_p, sizeof(int), 1, fd) ; // type_p
508  if (save_base) {
509  fwrite_be(colloc_r, sizeof(int), nzone, fd) ; // colloc_r
510  }
511  else
512  for (int l=0; l<nzone; l++)
513  if (colloc_r[l] != BASE_CHEB) {
514  cout << "Mg3d::sauve(FILE*, bool) : " << endl ;
515  cout << "The multi-grid is not with Chebyshev basis!!" << endl ;
516  cout << "Consider setting the 'save_base' flaf to 'true'!!"
517  << endl ;
518  arrete() ;
519 
520  }
521 }
522 
523  //--------------------------//
524  // Surcharge des operateurs //
525  //--------------------------//
526 
527 // Operateur <<
528 ostream& operator<<(ostream& o, const Mg3d& g) {
529  const char* tr[3] ;
530  tr[FIN] = "FIN" ; tr[RARE] = "RARE" ; tr[UNSURR] = "UNSURR" ;
531  const char* tang[2] ;
532  tang[NONSYM] = "NONSYM" ; tang[SYM] = "SYM" ;
533  const char* tbase[3] ;
534  tbase[BASE_CHEB] = "Chebyshev" ; tbase[BASE_LEG] = "Legendre" ;
535  tbase[BASE_JAC02] = "Jacobi(0,2)" ;
536  o << "Number of domains: " << g.nzone << endl ;
537  for (int i=0 ; i< g.nzone ; i++) {
538  o << " Domain #" << i << ": "
539  << "nr = " << g.nr[i] << ", " << tr[g.type_r[i]] << "; "
540  << "nt = " << g.nt[i] << ", " << tang[g.type_t] << "; "
541  << "np = " << g.np[i] << ", " << tang[g.type_p] << "; "
542  << "Collocation points type : " << tbase[g.colloc_r[i]] << endl ;
543  }
544  o << endl ;
545  return o ;
546 }
547 
548 // Operateur !=
549 bool Mg3d::operator!=(const Mg3d & titi) const {
550 
551  if (nzone != titi.nzone) return true ; // C'est vrai que c'est faux...
552 
553  for (int i=0 ; i<nzone ; i++) {
554  if (nr[i] != titi.nr[i]) return true ;
555  if (nt[i] != titi.nt[i]) return true ;
556  if (np[i] != titi.np[i]) return true ;
557 
558  if (type_r[i] != titi.type_r[i]) return true ;
559  if (colloc_r[i] != titi.colloc_r[i]) return true ;
560  }
561 
562  if (type_t != titi.type_t) return true ;
563  if (type_p != titi.type_p) return true ;
564 
565  // C'est faux que c'est vrai...
566  return false ;
567 }
568 
569 
570  //----------------------------------//
571  // Management of derived quantities //
572  //----------------------------------//
573 
574 void Mg3d::del_deriv() const {
575 
576  if (g_angu != 0x0) delete g_angu ;
577  if (g_angu_1dom != 0x0) delete g_angu_1dom ;
578  if (g_radial != 0x0) delete g_radial ;
579  if (g_twice != 0x0) delete g_twice ;
580  if (g_plus_half != 0x0) delete g_plus_half ;
581  if (g_plus_half_angu != 0x0) delete g_plus_half_angu ;
582  if (g_non_axi != 0x0) delete g_non_axi ;
583 
584  set_deriv_0x0() ;
585 
586 }
587 
588 void Mg3d::set_deriv_0x0() const {
589 
590  g_angu = 0x0 ;
591  g_angu_1dom = 0x0 ;
592  g_radial = 0x0 ;
593  g_twice = 0x0 ;
594  g_plus_half = 0x0 ;
595  g_plus_half_angu = 0x0 ;
596  g_non_axi = 0x0 ;
597 }
598 
599 
600  //--------------//
601  // Angular grid //
602  //--------------//
603 
604 const Mg3d* Mg3d::get_angu() const {
605 
606  if (g_angu == 0x0) { // The construction is required
607 
608  int* nbr_angu = new int[nzone] ;
609  for (int i=0 ; i<nzone ; i++) {
610  nbr_angu[i] = 1 ;
611  }
612  g_angu = new Mg3d(nzone, nbr_angu, type_r, nt, type_t, np, type_p,
613  colloc_r) ;
614  delete [] nbr_angu ;
615  }
616 
617  return g_angu ;
618 
619 }
620 
621  //-----------------------------//
622  // Angular grid for one domain //
623  //-----------------------------//
624 
625 const Mg3d* Mg3d::get_angu_1dom() const {
626 
627  if (g_angu_1dom == 0x0) { // The construction is required
628  int* nbr_angu = new int(1) ;
629  int* nbt_angu = new int(nt[0]) ;
630  int* nbp_angu = new int(np[0]) ;
631  int* type_r_angu = new int(FIN) ;
632 
633  g_angu_1dom = new Mg3d(1, nbr_angu, type_r_angu, nbt_angu, type_t,
634  nbp_angu, type_p) ;
635  delete nbr_angu ;
636  delete nbt_angu ;
637  delete nbp_angu ;
638  delete type_r_angu ;
639  }
640 
641  return g_angu_1dom ;
642 
643 }
644 
645  //--------------//
646  // Radial grid //
647  //--------------//
648 
649 const Mg3d* Mg3d::get_radial() const {
650 
651  if (g_radial == 0x0) { // The construction is required
652 
653  int* nbr_radial = new int[nzone] ;
654  for (int i=0 ; i<nzone ; i++) {
655  nbr_radial[i] = 1 ;
656  }
657  g_radial = new Mg3d(nzone, nr, type_r, nbr_radial, SYM, nbr_radial,
658  SYM, colloc_r) ;
659  delete [] nbr_radial ;
660  }
661 
662  return g_radial ;
663 
664 }
665 
666  //--------------------------------------//
667  // Grid with twice the number of points //
668  //--------------------------------------//
669 
670 const Mg3d* Mg3d::get_twice() const {
671 
672  if (g_twice == 0x0) { // The construction is required
673 
674  int* nbr = new int[nzone] ;
675  int* nbt = new int[nzone] ;
676  int* nbp = new int[nzone] ;
677 
678  for (int l=0; l<nzone; l++) {
679  if (nr[l] == 1) {
680  nbr[l] = 1 ;
681  }
682  else {
683  nbr[l] = 2*nr[l] - 1 ;
684  }
685 
686  if (nt[l] == 1) {
687  nbt[l] = 1 ;
688  }
689  else {
690  nbt[l] = 2*nt[l] - 1 ;
691  }
692 
693  if (np[l] == 1) {
694  nbp[l] = 1 ;
695  }
696  else {
697  nbp[l] = 2*np[l] ;
698  }
699  }
700 
701  g_twice = new Mg3d(nzone, nbr, type_r, nbt, type_t, nbp, type_p, colloc_r) ;
702 
703  delete [] nbr ;
704  delete [] nbt ;
705  delete [] nbp ;
706 
707  }
708 
709  return g_twice ;
710 
711 }
712 
713 
714  //--------------------------------------//
715  // Grid with 50% additional points in r //
716  //--------------------------------------//
717 
718 const Mg3d* Mg3d::plus_half() const {
719 
720  if (g_plus_half == 0x0) { // The construction is required
721 
722  int* nbr = new int[nzone] ;
723 
724  for (int l=0; l<nzone; l++) {
725  if (nr[l] == 1)
726  nbr[l] = 1 ;
727  else
728  nbr[l] = (3*nr[l])/2 ;
729  }
730 
731  g_plus_half = new Mg3d(nzone, nbr, type_r, nt, type_t, np, type_p, colloc_r) ;
732 
733  delete [] nbr ;
734 
735 
736  }
737 
738  return g_plus_half ;
739 
740 }
741 
742  //-------------------------------------------//
743  // Grid with 50% additional points in angles //
744  //-------------------------------------------//
745 
746 const Mg3d* Mg3d::plus_half_angu() const {
747 
748  if (g_plus_half_angu == 0x0) { // The construction is required
749 
750  int* nbt = new int[nzone] ;
751  int* nbp = new int[nzone] ;
752 
753  for (int l=0; l<nzone; l++) {
754  if (nt[l] == 1)
755  nbt[l] = 1 ;
756  else {
757  int nt_bis = (3*(nt[l]))/2;
758  if (nt_bis % 2 == 0) nt_bis += 1 ;
759  nbt[l] = nt_bis ;
760  }
761 
762  if (np[l] == 1)
763  nbp[l] = 1 ;
764  else
765  nbp[l] = 3*np[l]/2 ;
766  }
767 
768  g_plus_half_angu = new Mg3d(nzone, nr, type_r, nbt, type_t, nbp, type_p, colloc_r) ;
769 
770  delete [] nbt ;
771  delete [] nbp ;
772 
773 
774  }
775 
776  return g_plus_half_angu ;
777 
778 }
779 
780  //----------------------------------------------//
781  // Grid for rotations (5/4 points in theta/phi) //
782  //----------------------------------------------//
783 
784 const Mg3d* Mg3d::get_non_axi() const {
785 
786  if (g_non_axi == 0x0) { // The construction is required
787 
788  int* nbt = new int[nzone] ;
789  int* nbp = new int[nzone] ;
790 
791  for (int l=0; l<nzone; l++) {
792  if (nt[l] < 5)
793  nbt[l] = 5 ;
794  else
795  nbt[l] = nt[l] ;
796  if (np[l] < 4)
797  nbp[l] = 4 ;
798  else
799  nbp[l] = np[l] ;
800  }
801 
802  g_non_axi = new Mg3d(nzone, nr, type_r, nbt, type_t, nbp, type_p, colloc_r) ;
803 
804  delete [] nbt ;
805  delete [] nbp ;
806 
807 
808  }
809 
810  return g_non_axi ;
811 
812 }
813 
814 
815 bool Mg3d::operator==(const Mg3d& mgi) const {
816 
817  bool resu = true ;
818 
819  if (mgi.get_nzone() != nzone) {
820  resu = false ;
821  }
822  else {
823  for (int i=0; i<nzone; i++) {
824  if (mgi.get_nr(i) != nr[i]) resu = false ;
825  if (mgi.get_np(i) != np[i]) resu = false ;
826  if (mgi.get_nt(i) != nt[i]) resu = false ;
827  if (mgi.get_type_r(i) != type_r[i]) resu =false ;
828  if (mgi.get_colloc_r(i) != colloc_r[i]) resu = false ;
829  }
830  }
831 
832  if (mgi.get_type_t() != type_t) resu = false ;
833  if (mgi.get_type_p() != type_p) resu = false ;
834 
835  return resu ;
836 
837 }
838 }
void set_deriv_0x0() const
Sets to 0x0 all the pointers on derived quantities (g_radial , g_angu, g_twice, ...
Definition: mg3d.C:588
int get_type_p() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the transformatio...
Definition: grilles.h:512
const Mg3d * get_non_axi() const
Returns the pointer on the grid which has at least 4 points in the direction and at least 5 in the ...
Definition: mg3d.C:784
Mg3d * g_non_axi
Pointer on the grid which has at least 4 points in the direction and at least 5 in the direction (f...
Definition: grilles.h:334
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:479
int type_p
Type of sampling in (SYM, NONSYM)
Definition: grilles.h:299
Lorene prototypes.
Definition: app_hor.h:67
Mg3d * g_plus_half
Pointer on the grid which has 50% more points in r dimension (for desaliasing).
Definition: grilles.h:323
int * np
Array (size: nzone) of nb. of points in .
Definition: grilles.h:288
Mg3d * g_angu_1dom
Pointer on the associated angular grid with only one domain.
Definition: grilles.h:312
const Mg3d * get_twice() const
Returns the pointer on the grid which has twice the number of points in each dimension (for desaliasi...
Definition: mg3d.C:670
void del_deriv() const
Deletes all the derived quantities (g_radial , g_angu, g_twice, ...)
Definition: mg3d.C:574
int get_type_t() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the equatorial pl...
Definition: grilles.h:502
const Mg3d * get_radial() const
Returns the pointer on the associated radial grid.
Definition: mg3d.C:649
bool operator==(const Mg3d &) const
Comparison operator (egality)
Definition: mg3d.C:815
int type_t
Type of sampling in (SYM, NONSYM)
Definition: grilles.h:296
int get_colloc_r(int l) const
Returns the type of collocation points used in domain no.
Definition: grilles.h:528
~Mg3d()
Destructor.
Definition: mg3d.C:480
Mg3d * g_plus_half_angu
Pointer on the grid which has 50% more points in and dimension (for desaliasing).
Definition: grilles.h:328
void sauve(FILE *fd, bool save_base=false) const
Saves into a file.
Definition: mg3d.C:500
const Mg3d * get_angu() const
Returns the pointer on the associated angular grid.
Definition: mg3d.C:604
bool search_file(ifstream &infile, const string &pattern)
A function that searches for a pattern in a file and places the file stream after the found pattern...
Definition: misc.C:53
int nzone
Number of domains (zones)
Definition: grilles.h:284
Mg3d * g_angu
Pointer on the associated angular grid.
Definition: grilles.h:310
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
int fwrite_be(const int *aa, int size, int nb, FILE *fich)
Writes integer(s) into a binary file according to the big endian convention.
Definition: fwrite_be.C:73
Mg3d * g_twice
Pointer on the grid which has twice the number of points in each dimension (for desaliasing).
Definition: grilles.h:318
Mg3d * g_radial
Pointer on the associated radial grid.
Definition: grilles.h:313
const Mg3d * plus_half() const
Returns the pointer on the grid which has 50% more points in r dimension (for desaliasing).
Definition: mg3d.C:718
Grille3d ** g
Array (size: nzone) of pointers on the Grille3d&#39;s.
Definition: grilles.h:308
int * type_r
Array (size: nzone) of type of sampling in r ( ) (RARE,FIN, UNSURR)
Definition: grilles.h:293
int fread_be(int *aa, int size, int nb, FILE *fich)
Reads integer(s) from a binary file according to the big endian convention.
Definition: fread_be.C:72
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:469
int * colloc_r
Array (size: nzone) of type of collocation points in r ( ) and related decompoisition bases (BASE_CHE...
Definition: grilles.h:305
Multi-domain grid.
Definition: grilles.h:279
const Mg3d * plus_half_angu() const
Returns the pointer on the grid which has 50% more points in and dimension (for desaliasing)...
Definition: mg3d.C:746
void arrete(int a=0)
Setting a stop point in a code.
Definition: arrete.C:64
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:474
bool operator!=(const Mg3d &) const
Operator !=.
Definition: mg3d.C:549
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition: grilles.h:491
Mg3d(int nz, int nbr[], int typr[], int nbt[], int typt, int nbp[], int typp, int *base_r=0x0)
General constructor.
Definition: mg3d.C:172
int * nt
Array (size: nzone) of nb. of points in .
Definition: grilles.h:287
const Mg3d * get_angu_1dom() const
Returns the pointer on the associated mono-domain angular grid.
Definition: mg3d.C:625
3D grid class in one domain.
Definition: grilles.h:200
int * nr
Array (size: nzone) of nb. of points in r ( )
Definition: grilles.h:286