LORENE
diff_x3dsdx.C
1/*
2 * Methods for the class Diff_x3dsdx
3 *
4 * (see file diff.h for documentation).
5 *
6 */
7
8/*
9 * Copyright (c) 2005 Jerome Novak
10 *
11 * This file is part of LORENE.
12 *
13 * LORENE is free software; you can redistribute it and/or modify
14 * it under the terms of the GNU General Public License version 2
15 * as published by the Free Software Foundation.
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
28char diff_x3dsdx_C[] = "$Header: /cvsroot/Lorene/C++/Source/Diff/diff_x3dsdx.C,v 1.5 2014/10/13 08:52:51 j_novak Exp $" ;
29
30/*
31 * $Id: diff_x3dsdx.C,v 1.5 2014/10/13 08:52:51 j_novak Exp $
32 * $Log: diff_x3dsdx.C,v $
33 * Revision 1.5 2014/10/13 08:52:51 j_novak
34 * Lorene classes and functions now belong to the namespace Lorene.
35 *
36 * Revision 1.4 2014/10/06 15:13:05 j_novak
37 * Modified #include directives to use c++ syntax.
38 *
39 * Revision 1.3 2007/12/11 15:28:11 jl_cornou
40 * Jacobi(0,2) polynomials partially implemented
41 *
42 * Revision 1.2 2005/01/11 15:52:38 j_novak
43 * Corrected a mistake...
44 *
45 * Revision 1.1 2005/01/11 15:16:10 j_novak
46 * More Diff operators.
47 *
48 * Revision 1.1 2005/01/10 16:34:52 j_novak
49 * New class for 1D mono-domain differential operators.
50 *
51 *
52 * $Header: /cvsroot/Lorene/C++/Source/Diff/diff_x3dsdx.C,v 1.5 2014/10/13 08:52:51 j_novak Exp $
53 *
54 */
55
56// C headers
57#include <cassert>
58#include <cstdlib>
59
60// Lorene headers
61#include "diff.h"
62#include "proto.h"
63
64namespace Lorene {
65void xpundsdx_1d(int, double**, int) ;
66void mult2_xp1_1d(int, double**, int) ;
67
68namespace {
69 int nap = 0 ;
70 Matrice* tab[MAX_BASE*Diff::max_points] ;
71 int nr_done[Diff::max_points] ;
72}
73
74Diff_x3dsdx::Diff_x3dsdx(int base_r, int nr) : Diff(base_r, nr) {
75 initialize() ;
76}
77
78Diff_x3dsdx::Diff_x3dsdx(const Diff_x3dsdx& diff_in) : Diff(diff_in) {
79 assert (nap != 0) ;
80}
81
83
85 if (nap == 0) {
86 for (int i=0; i<max_points; i++) {
87 nr_done[i] = -1 ;
88 for (int j=0; j<MAX_BASE; j++)
89 tab[j*max_points+i] = 0x0 ;
90 }
91 nap = 1 ;
92 }
93 return ;
94}
95
96void Diff_x3dsdx::operator=(const Diff_x3dsdx& diff_in) {
97 assert (nap != 0) ;
98 Diff::operator=(diff_in) ;
99
100}
101
103
104 bool done = false ;
105 int indice ;
106 for (indice =0; indice<max_points; indice++) {
107 if (nr_done[indice] == npoints) {
108 if (tab[base*max_points + indice] != 0x0) done = true ;
109 break ;
110 }
111 if (nr_done[indice] == -1)
112 break ;
113 }
114 if (!done) { //The computation must be done ...
115 if (indice == max_points) {
116 cerr << "Diff_x3dsdx::get_matrice() : no space left!!" << '\n'
117 << "The value of Diff.max_points must be increased..." << endl ;
118 abort() ;
119 }
120 nr_done[indice] = npoints ;
121 tab[base*max_points + indice] = new Matrice(npoints, npoints) ;
122 Matrice& resu = *tab[base*max_points + indice] ;
123 resu.set_etat_qcq() ;
124
125 double* vect = new double[npoints] ;
126 double* cres = new double[npoints] ;
127 for (int i=0; i<npoints; i++) {
128 for (int j=0; j<npoints; j++)
129 vect[j] = 0. ;
130 vect[i] = 1. ;
131 if (base == R_CHEBU) {
132 dsdx_1d(npoints, &vect, R_CHEBU) ;
133 mult_xm1_1d_cheb(npoints, vect, cres) ;
134 mult2_xm1_1d_cheb(npoints, cres, vect) ;
135 for (int j=0; j<npoints; j++)
136 resu.set(j,i) = vect[j] ;
137 }
138 else if (base == R_JACO02) {
139 xpundsdx_1d(npoints, &vect, base << TRA_R) ;
140 mult2_xp1_1d(npoints, &vect, base << TRA_R) ;
141 for (int j=0; j<npoints ; j++)
142 resu.set(j,i) = vect[j] ;
143 }
144 else {
145 xdsdx_1d(npoints, &vect, base << TRA_R) ;
146 multx2_1d(npoints, &vect, base << TRA_R) ;
147 for (int j=0; j<npoints; j++)
148 resu.set(j,i) = vect[j] ;
149
150 }
151 }
152 delete [] vect ;
153 delete [] cres ;
154 }
155
156 return *tab[base*max_points + indice] ;
157}
158
159ostream& Diff_x3dsdx::operator>>(ostream& ost) const {
160
161 ost << " xi^3 * d / dx " << endl ;
162 return ost ;
163
164}
165}
Class for the elementary differential operator (see the base class Diff ).
Definition diff.h:611
void initialize()
Initializes arrays.
Definition diff_x3dsdx.C:84
virtual ostream & operator>>(ostream &) const
Operator >> (virtual function called by the operator <<).
Diff_x3dsdx(int base_r, int nr)
Standard constructor.
Definition diff_x3dsdx.C:74
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
void operator=(const Diff_x3dsdx &)
Assignment to another Diff_x3dsdx.
Definition diff_x3dsdx.C:96
virtual ~Diff_x3dsdx()
Destructor.
Definition diff_x3dsdx.C:82
Base (abstract) class for 1D spectral differential operators in one domain.
Definition diff.h:65
int npoints
Number of coefficients.
Definition diff.h:75
void operator=(const Diff &)
Assignment to another Diff.
Definition diff.C:75
static const int max_points
Maximal number of matrices stored per base.
Definition diff.h:71
int base
Base in radial direction.
Definition diff.h:74
Matrix handling.
Definition matrice.h:152
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition matrice.C:175
double & set(int j, int i)
Read/write of a particuliar element.
Definition matrice.h:277
#define MAX_BASE
Nombre max. de bases differentes.
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
Lorene prototypes.
Definition app_hor.h:64