ParMooN
 All Classes Functions Variables Friends Pages
BSExample.h
1 // Navier-Stokes problem, solution in ansatz space
2 // velocity pw quadratic, pressure linear
3 //
4 //
5 void ExampleFile()
6 {
7  OutPut("Example: BSExample.h" << endl) ;
8 }
9 
10 // ========================================================================
11 // exact solution
12 // ========================================================================
13 void ExactU1(double x, double y, double z, double *values)
14 {
15  values[0] = sin(Pi*x)*sin(Pi*y)*sin(Pi*z)+x*x*x*x*cos(Pi*y);
16  values[1] = cos(Pi*x)*Pi*sin(Pi*y)*sin(Pi*z)+4*x*x*x*cos(Pi*y);
17  values[2] = sin(Pi*x)*cos(Pi*y)*Pi*sin(Pi*z)-x*x*x*x*sin(Pi*y)*Pi;
18  values[3] = sin(Pi*x)*sin(Pi*y)*cos(Pi*z)*Pi;
19  values[4] = -3*sin(Pi*x)*Pi*Pi*sin(Pi*y)*sin(Pi*z)
20  +12*x*x*cos(Pi*y)-x*x*x*x*cos(Pi*y)*Pi*Pi;
21 }
22 
23 void ExactU2(double x, double y, double z, double *values)
24 {
25  values[0] = cos(Pi*x)*cos(Pi*y)*cos(Pi*z)-3*y*y*y*z;
26  values[1] = -Pi*sin(Pi*x)*cos(Pi*y)*cos(Pi*z);
27  values[2] = -Pi*cos(Pi*x)*sin(Pi*y)*cos(Pi*z)-9*y*y*z;
28  values[3] = -Pi*cos(Pi*x)*cos(Pi*y)*sin(Pi*z)-3*y*y*y;
29  values[4] = -3*cos(Pi*x)*Pi*Pi*cos(Pi*y)*cos(Pi*z)-18*y*z;
30 }
31 
32 void ExactU3(double x, double y, double z, double *values)
33 {
34  values[0] = cos(Pi*x)*sin(Pi*y)*cos(Pi*z)+cos(Pi*x)*sin(Pi*y)*sin(Pi*z)
35  -4*x*x*x*cos(Pi*y)*z+4.5*y*y*z*z;
36  values[1] = -sin(Pi*x)*sin(Pi*y)*cos(Pi*z)*Pi
37  -12*x*x*cos(Pi*y)*z-sin(Pi*z)*sin(Pi*x)*Pi*sin(Pi*y);
38  values[2] = cos(Pi*z)*cos(Pi*x)*cos(Pi*y)*Pi
39  +4*x*x*x*sin(Pi*y)*Pi*z+cos(Pi*x)*cos(Pi*y)*sin(Pi*z)*Pi+9*y*z*z;
40  values[3] = -cos(Pi*x)*Pi*sin(Pi*y)*sin(Pi*z)
41  -4*x*x*x*cos(Pi*y)+cos(Pi*x)*sin(Pi*y)*Pi*cos(Pi*z)+9*y*y*z;
42  values[4] = -3*cos(Pi*x)*Pi*Pi*sin(Pi*y)*cos(Pi*z)
43  -24*x*cos(Pi*y)*z
44  -3*sin(Pi*z)*cos(Pi*x)*Pi*Pi*sin(Pi*y)
45  +4*x*x*x*cos(Pi*y)*Pi*Pi*z
46  +9*z*z+9*y*y;
47 }
48 
49 void ExactP(double x, double y, double z, double *values)
50 {
51  values[0] = 3*x-sin(y+4*z)-1.5-sin(1.0)*cos(1.0)
52  -4*cos(1.0)*cos(1.0)*cos(1.0)*cos(1.0)*sin(1.0)
53  +sin(1.0)*cos(1.0)*cos(1.0)-2*sin(1.0)*sin(1.0)*sin(1.0)
54  +2*sin(1.0)*cos(1.0)*cos(1.0)*cos(1.0)+2*sin(1.0);
55  values[1] = 3;
56  values[2] = -cos(y+4*z);
57  values[3] = -4*cos(y+4*z);
58  values[4] = 0;
59 }
60 
61 // kind of boundary condition (for FE space needed)
62 void BoundCondition(int CompID, double x, double y, double z, BoundCond &cond)
63 {
64  cond = DIRICHLET;
65 // cond = NEUMANN;
66 }
67 
68 // value of boundary condition
69 void U1BoundValue(int CompID, double x, double y, double z, double &value)
70 {
71  value = sin(Pi*x)*sin(Pi*y)*sin(Pi*z)+x*x*x*x*cos(Pi*y);
72 }
73 
74 // value of boundary condition
75 void U2BoundValue(int CompID, double x, double y, double z, double &value)
76 {
77  value = cos(Pi*x)*cos(Pi*y)*cos(Pi*z)-3*y*y*y*z;
78 }
79 
80 // value of boundary condition
81 void U3BoundValue(int CompID, double x, double y, double z, double &value)
82 {
83  value = cos(Pi*x)*sin(Pi*y)*cos(Pi*z)+cos(Pi*x)*sin(Pi*y)*sin(Pi*z)
84  -4*x*x*x*cos(Pi*y)*z+4.5*y*y*z*z;
85 }
86 
87 // ========================================================================
88 // exact solution
89 // ========================================================================
90 void InitialU1(double x, double y, double z, double *values)
91 {
92  values[0] = sin(Pi*x)*sin(Pi*y)*sin(Pi*z)+x*x*x*x*cos(Pi*y);
93  values[1] = cos(Pi*x)*Pi*sin(Pi*y)*sin(Pi*z)+4*x*x*x*cos(Pi*y);
94  values[2] = sin(Pi*x)*cos(Pi*y)*Pi*sin(Pi*z)-x*x*x*x*sin(Pi*y)*Pi;
95  values[3] = sin(Pi*x)*sin(Pi*y)*cos(Pi*z)*Pi;
96  values[4] = -3*sin(Pi*x)*Pi*Pi*sin(Pi*y)*sin(Pi*z)
97  +12*x*x*cos(Pi*y)-x*x*x*x*cos(Pi*y)*Pi*Pi;
98  values[0] = 0;
99 }
100 
101 void InitialU2(double x, double y, double z, double *values)
102 {
103  values[0] = cos(Pi*x)*cos(Pi*y)*cos(Pi*z)-3*y*y*y*z;
104  values[1] = -Pi*sin(Pi*x)*cos(Pi*y)*cos(Pi*z);
105  values[2] = -Pi*cos(Pi*x)*sin(Pi*y)*cos(Pi*z)-9*y*y*z;
106  values[3] = -Pi*cos(Pi*x)*cos(Pi*y)*sin(Pi*z)-3*y*y*y;
107  values[4] = -3*cos(Pi*x)*Pi*Pi*cos(Pi*y)*cos(Pi*z)-18*y*z;
108  values[0] = 0;
109 }
110 
111 void InitialU3(double x, double y, double z, double *values)
112 {
113  values[0] = cos(Pi*x)*sin(Pi*y)*cos(Pi*z)+cos(Pi*x)*sin(Pi*y)*sin(Pi*z)
114  -4*x*x*x*cos(Pi*y)*z+4.5*y*y*z*z;
115  values[1] = -sin(Pi*x)*sin(Pi*y)*cos(Pi*z)*Pi
116  -12*x*x*cos(Pi*y)*z-sin(Pi*z)*sin(Pi*x)*Pi*sin(Pi*y);
117  values[2] = cos(Pi*z)*cos(Pi*x)*cos(Pi*y)*Pi
118  +4*x*x*x*sin(Pi*y)*Pi*z+cos(Pi*x)*cos(Pi*y)*sin(Pi*z)*Pi+9*y*z*z;
119  values[3] = -cos(Pi*x)*Pi*sin(Pi*y)*sin(Pi*z)
120  -4*x*x*x*cos(Pi*y)+cos(Pi*x)*sin(Pi*y)*Pi*cos(Pi*z)+9*y*y*z;
121  values[4] = -3*cos(Pi*x)*Pi*Pi*sin(Pi*y)*cos(Pi*z)
122  -24*x*cos(Pi*y)*z
123  -3*sin(Pi*z)*cos(Pi*x)*Pi*Pi*sin(Pi*y)
124  +4*x*x*x*cos(Pi*y)*Pi*Pi*z
125  +9*z*z+9*y*y;
126  values[0] = 0;
127 }
128 
129 void InitialP(double x, double y, double z, double *values)
130 {
131  values[0] = 3*x-sin(y+4*z)-1.5-sin(1.0)*cos(1.0)
132  -4*cos(1.0)*cos(1.0)*cos(1.0)*cos(1.0)*sin(1.0)
133  +sin(1.0)*cos(1.0)*cos(1.0)-2*sin(1.0)*sin(1.0)*sin(1.0)
134  +2*sin(1.0)*cos(1.0)*cos(1.0)*cos(1.0)+2*sin(1.0);
135  values[1] = 3;
136  values[2] = -cos(y+4*z);
137  values[3] = -4*cos(y+4*z);
138  values[4] = 0;
139  values[0] = 0;
140 }
141 
142 // ========================================================================
143 // coefficients for Stokes form: A, B1, B2, f1, f2
144 // ========================================================================
145 void LinCoeffs(int n_points, double *X, double *Y, double *Z,
146  double **parameters, double **coeffs)
147 {
148  double eps = 1./TDatabase::ParamDB->RE_NR;
149  int i;
150  double *coeff, x, y, z, u1, u2, u3, ux, uy, uz;
151 
152  if (TDatabase::ParamDB->FLOW_PROBLEM_TYPE==STOKES)
153  {
154  for(i=0;i<n_points;i++)
155  {
156  coeff = coeffs[i];
157 
158  x = X[i];
159  y = Y[i];
160  z = Z[i];
161  coeff[0] = eps;
162  coeff[1] = -eps*(-3*sin(Pi*x)*Pi*Pi*sin(Pi*y)*sin(Pi*z)
163  +12*x*x*cos(Pi*y)-x*x*x*x*cos(Pi*y)*Pi*Pi) + 3; // f1
164  coeff[2] = -eps*(-3*cos(Pi*x)*Pi*Pi*cos(Pi*y)*cos(Pi*z)-18*y*z)
165  - cos(y+4*z); // f2
166  coeff[3] = -eps*(-3*cos(Pi*x)*Pi*Pi*sin(Pi*y)*cos(Pi*z)
167  -24*x*cos(Pi*y)*z
168  -3*sin(Pi*z)*cos(Pi*x)*Pi*Pi*sin(Pi*y)
169  +4*x*x*x*cos(Pi*y)*Pi*Pi*z
170  +9*z*z+9*y*y)
171  - 4*cos(y+4*z); // f3
172  }
173  }
174  else if(TDatabase::ParamDB->FLOW_PROBLEM_TYPE==NSE)
175  {
176  for(i=0;i<n_points;i++)
177  {
178  coeff = coeffs[i];
179 
180  x = X[i];
181  y = Y[i];
182  z = Z[i];
183  u1 = sin(Pi*x)*sin(Pi*y)*sin(Pi*z)+x*x*x*x*cos(Pi*y);
184  u2 = cos(Pi*x)*cos(Pi*y)*cos(Pi*z)-3*y*y*y*z;
185  u3 = cos(Pi*x)*sin(Pi*y)*cos(Pi*z)+cos(Pi*x)*sin(Pi*y)*sin(Pi*z)
186  -4*x*x*x*cos(Pi*y)*z+4.5*y*y*z*z;
187  ux = cos(Pi*x)*Pi*sin(Pi*y)*sin(Pi*z)+4*x*x*x*cos(Pi*y);
188  uy = sin(Pi*x)*cos(Pi*y)*Pi*sin(Pi*z)-x*x*x*x*sin(Pi*y)*Pi;
189  uz = sin(Pi*x)*sin(Pi*y)*cos(Pi*z)*Pi;
190  coeff[1] = -eps*(-3*sin(Pi*x)*Pi*Pi*sin(Pi*y)*sin(Pi*z)
191  +12*x*x*cos(Pi*y)-x*x*x*x*cos(Pi*y)*Pi*Pi)
192  + u1*ux+u2*uy+u3*uz+ 3; // f1
193 
194  ux = -Pi*sin(Pi*x)*cos(Pi*y)*cos(Pi*z);
195  uy = -Pi*cos(Pi*x)*sin(Pi*y)*cos(Pi*z)-9*y*y*z;
196  uz= -Pi*cos(Pi*x)*cos(Pi*y)*sin(Pi*z)-3*y*y*y;
197  coeff[2] = -eps*(-3*cos(Pi*x)*Pi*Pi*cos(Pi*y)*cos(Pi*z)-18*y*z)
198  + u1*ux+u2*uy+u3*uz - cos(y+4*z); // f2
199 
200  ux = -sin(Pi*x)*sin(Pi*y)*cos(Pi*z)*Pi
201  -12*x*x*cos(Pi*y)*z-sin(Pi*z)*sin(Pi*x)*Pi*sin(Pi*y);
202  uy = cos(Pi*z)*cos(Pi*x)*cos(Pi*y)*Pi
203  +4*x*x*x*sin(Pi*y)*Pi*z+cos(Pi*x)*cos(Pi*y)*sin(Pi*z)*Pi+9*y*z*z;
204  uz = -cos(Pi*x)*Pi*sin(Pi*y)*sin(Pi*z)
205  -4*x*x*x*cos(Pi*y)+cos(Pi*x)*sin(Pi*y)*Pi*cos(Pi*z)+9*y*y*z;
206  coeff[3] = -eps*(-3*cos(Pi*x)*Pi*Pi*sin(Pi*y)*cos(Pi*z)
207  -24*x*cos(Pi*y)*z
208  -3*sin(Pi*z)*cos(Pi*x)*Pi*Pi*sin(Pi*y)
209  +4*x*x*x*cos(Pi*y)*Pi*Pi*z
210  +9*z*z+9*y*y)
211  + u1*ux+u2*uy+u3*uz - 4*cos(y+4*z); // f3
212  coeff[0] = eps;
213  }
214  }
215  else
216  { OutPut("Oseen type not included.... Change FLOW_PROBLEM_TYPE"<<endl);
217  exit(0);
218 
219  }
220 }
TParaDB::RE_NR
double RE_NR
Definition: Database.h:313
TDatabase::ParamDB
static TParamDB * ParamDB
Definition: Database.h:1134
Generated on Tue May 10 2016 09:07:56 for ParMooN by   doxygen 1.8.5