OpenLB 1.7
Loading...
Searching...
No Matches
frameChangeF3D.hh
Go to the documentation of this file.
1/* This file is part of the OpenLB library
2 *
3 * Copyright (C) 2013, 2015 Gilles Zahnd, Mathias J. Krause
4 * Marie-Luise Maier
5 * E-mail contact: info@openlb.net
6 * The most recent release of OpenLB can be downloaded at
7 * <http://www.openlb.net/>
8 *
9 * This program is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU General Public License
11 * as published by the Free Software Foundation; either version 2
12 * of the License, or (at your option) any later version.
13 *
14 * This program is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 * GNU General Public License for more details.
18 *
19 * You should have received a copy of the GNU General Public
20 * License along with this program; if not, write to the Free
21 * Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
22 * Boston, MA 02110-1301, USA.
23*/
24
25#ifndef FRAME_CHANGE_F_3D_HH
26#define FRAME_CHANGE_F_3D_HH
27
28#include<cmath>
29#include <random>
30
31#include "frameChangeF3D.h"
32#include "frameChangeF2D.h"
33#include "core/superLattice.h"
34#include "dynamics/lbm.h" // for computation of lattice rho and velocity
35#include "utilities/vectorHelpers.h" // for normalize
36#include "geometry/superGeometry.h"
37
38#ifndef M_PI
39#define M_PI 3.14159265358979323846
40#endif
41
42namespace olb {
43
44
45template <typename T>
46RotatingLinear3D<T>::RotatingLinear3D(std::vector<T> axisPoint_,
47 std::vector<T> axisDirection_, T w_, T scale_)
48 : AnalyticalF3D<T,T>(3), axisPoint(axisPoint_), axisDirection(axisDirection_),
49 w(w_), scale(scale_) { }
50
51
52template <typename T>
53bool RotatingLinear3D<T>::operator()(T output[], const T x[])
54{
55 output[0] = (axisDirection[1]*(x[2]-axisPoint[2]) - axisDirection[2]*(x[1]-axisPoint[1]))*w*scale;
56 output[1] = (axisDirection[2]*(x[0]-axisPoint[0]) - axisDirection[0]*(x[2]-axisPoint[2]))*w*scale;
57 output[2] = (axisDirection[0]*(x[1]-axisPoint[1]) - axisDirection[1]*(x[0]-axisPoint[0]))*w*scale;
58 return true;
59}
60
61
62
63template <typename T>
64RotatingLinearAnnulus3D<T>::RotatingLinearAnnulus3D(std::vector<T> axisPoint_,
65 std::vector<T> axisDirection_, T w_, T ri_, T ro_, T scale_)
66 : AnalyticalF3D<T,T>(3), axisPoint(axisPoint_), axisDirection(axisDirection_),
67 w(w_), ri(ri_), ro(ro_), scale(scale_) { }
68
69
70template <typename T>
71bool RotatingLinearAnnulus3D<T>::operator()(T output[], const T x[])
72{
73
74 if ( (util::sqrt((x[0]-axisPoint[0])*(x[0]-axisPoint[0])+(x[1]-axisPoint[1])*(x[1]-axisPoint[1])) < ri) || (util::sqrt((x[0]-axisPoint[0])*(x[0]-axisPoint[0])+(x[1]-axisPoint[1])*(x[1]-axisPoint[1])) >= ro) ) {
75 output[0] = 0.;
76 output[1] = 0.;
77 output[2] = 0.;
78 }
79 else {
80 T L[3];
81 T si[3];
82 T so[3];
83 int sign1 = 1;
84 int sign2 = 1;
85 if ( (axisDirection[2] && (x[0] == axisPoint[0])) || (axisDirection[1] && (x[2] == axisPoint[2])) || (axisDirection[0] && (x[1] == axisPoint[1])) ) {
86 int sign = 1;
87 if ( (axisDirection[2] && (x[1] < axisPoint[1])) || (axisDirection[1] && (x[0] < axisPoint[0])) || (axisDirection[0] && (x[2] < axisPoint[2])) ) {
88 sign = -1;
89 }
90 //Compute point of intersection between the inner cylinder and the line between axisPoint and x
91 si[0] = axisDirection[0]*x[0] + axisDirection[1]*(axisPoint[0]+sign*ri) + axisDirection[2]*x[0];
92 si[1] = axisDirection[0]*x[1] + axisDirection[1]*x[1] + axisDirection[2]*(axisPoint[1]+sign*ri);
93 si[2] = axisDirection[0]*(axisPoint[2]+sign*ri) + axisDirection[1]*x[2] + axisDirection[2]*x[2];
94 //Compute point of intersection between the outer cylinder and the line between axisPoint and x
95 so[0] = axisDirection[0]*x[0] + axisDirection[1]*(axisPoint[0]+sign*ro) + axisDirection[2]*x[0];
96 so[1] = axisDirection[0]*x[1] + axisDirection[1]*x[1] + axisDirection[2]*(axisPoint[1]+sign*ro);
97 so[2] = axisDirection[0]*(axisPoint[2]+sign*ro) + axisDirection[1]*x[2] + axisDirection[2]*x[2];
98 }
99 else {
100 T alpha;
101 //which quadrant
102 if ( (axisDirection[2] && (x[0] < axisPoint[0])) || (axisDirection[1] && (x[2] < axisPoint[2])) || (axisDirection[0] && (x[1] < axisPoint[1])) ) {
103 sign1 = -1;
104 }
105 if ( (axisDirection[2] && (x[1] < axisPoint[1])) || (axisDirection[1] && (x[0] < axisPoint[0])) || (axisDirection[0] && (x[2] < axisPoint[2])) ) {
106 sign2 = -1;
107 }
108 alpha = util::atan( ( sign2*(axisDirection[0]*(x[2]-axisPoint[2]) + axisDirection[1]*(x[0]-axisPoint[0]) + axisDirection[2]*(x[1]-axisPoint[1]) ) ) / \
109 (sign1*(axisDirection[0]*(x[1]-axisPoint[1]) + axisDirection[1]*(x[2]-axisPoint[2]) + axisDirection[2]*(x[0]-axisPoint[0]))) );
110 si[0] = axisDirection[0]*x[0] + axisDirection[1]*sign2*util::sin(alpha)*ri + axisDirection[2]*sign1*util::cos(alpha)*ri;
111 si[1] = axisDirection[0]*sign1*util::cos(alpha)*ri + axisDirection[1]*x[1] + axisDirection[2]*sign2*util::sin(alpha)*ri;
112 si[2] = axisDirection[0]*sign2*util::sin(alpha)*ri + axisDirection[1]*sign1*util::cos(alpha)*ri + axisDirection[2]*x[2];
113 so[0] = axisDirection[0]*x[0] + axisDirection[1]*sign2*util::sin(alpha)*ro + axisDirection[2]*sign1*util::cos(alpha)*ro;
114 so[1] = axisDirection[0]*sign1*util::cos(alpha)*ro + axisDirection[1]*x[1] + axisDirection[2]*sign2*util::sin(alpha)*ro;
115 so[2] = axisDirection[0]*sign2*util::sin(alpha)*ro + axisDirection[1]*sign1*util::cos(alpha)*ro + axisDirection[2]*x[2];
116 }
117
118 //Compute difference of intersections in all directions
119 L[0] = so[0]-si[0];
120 L[1] = so[1]-si[1];
121 L[2] = so[2]-si[2];
122 bool b0 = (L[0] == 0.);
123 bool b1 = (L[1] == 0.);
124 bool b2 = (L[2] == 0.);
125
126 output[0] = ((axisDirection[1]*(axisPoint[2]-si[2]) - axisDirection[2]*(axisPoint[1]-si[1])) *
127 (1 - (axisDirection[1]*(x[2]-si[2])/(L[2]+b2) + axisDirection[2]*(x[1]-si[1])/(L[1]+b1))) )*w*scale;
128 output[1] = ((axisDirection[2]*(axisPoint[0]-si[0]) - axisDirection[0]*(axisPoint[2]-si[2])) *
129 (1 - (axisDirection[2]*(x[0]-si[0])/(L[0]+b0) + axisDirection[0]*(x[2]-si[2])/(L[2]+b2))) )*w*scale;
130 output[2] = ((axisDirection[0]*(axisPoint[1]-si[1]) - axisDirection[1]*(axisPoint[0]-si[0])) *
131 (1 - (axisDirection[0]*(x[1]-si[1])/(L[1]+b1) + axisDirection[1]*(x[0]-si[0])/(L[0]+b0))) )*w*scale;
132
133 }
134 return true;
135}
136
137template <typename T>
138RotatingQuadratic1D<T>::RotatingQuadratic1D(std::vector<T> axisPoint_,
139 std::vector<T> axisDirection_, T w_, T scale_, T additive_)
140 : AnalyticalF3D<T,T>(1), w(w_), scale(scale_), additive(additive_)
141{
142 axisPoint.resize(3);
143 axisDirection.resize(3);
144 for (int i = 0; i < 3; ++i) {
145 axisPoint[i] = axisPoint_[i];
146 axisDirection[i] = axisDirection_[i];
147 }
148}
149
150
151template <typename T>
152bool RotatingQuadratic1D<T>::operator()(T output[], const T x[])
153{
154 T radiusSquared = ((x[1]-axisPoint[1])*(x[1]-axisPoint[1])
155 +(x[2]-axisPoint[2])*(x[2]-axisPoint[2]))*axisDirection[0]
156 + ((x[0]-axisPoint[0])*(x[0]-axisPoint[0])
157 +(x[2]-axisPoint[2])*(x[2]-axisPoint[2]))*axisDirection[1]
158 + ((x[1]-axisPoint[1])*(x[1]-axisPoint[1])
159 +(x[0]-axisPoint[0])*(x[0]-axisPoint[0]))*axisDirection[2];
160 output[0] = scale*w*w/2*radiusSquared+additive;
161 return true;
162}
163
164
165template <typename T>
166CirclePowerLaw3D<T>::CirclePowerLaw3D(olb::Vector<T, 3> center, std::vector<T> normal,
167 T maxVelocity, T radius, T n, T scale)
168 : AnalyticalF3D<T,T>(3), _center(center), _normal(util::normalize(normal)),
169 _radius(radius), _maxVelocity(maxVelocity), _n(n), _scale(scale) { }
170
171template <typename T>
172CirclePowerLaw3D<T>::CirclePowerLaw3D(T center0, T center1, T center2, T normal0, T normal1, T normal2, T radius, T maxVelocity, T n, T scale )
173 : AnalyticalF3D<T,T>(3), _radius(radius), _maxVelocity(maxVelocity), _n(n), _scale(scale)
174{
175// _center.push_back(center0);
176 _center[0] = center0;
177 _center[1] = center1;
178 _center[2] = center2;
179 std::vector<T> normalTmp;
180 normalTmp.push_back(normal0);
181 normalTmp.push_back(normal1);
182 normalTmp.push_back(normal2 );
183 _normal = normalTmp;
184}
185
186template <typename T>
187CirclePowerLaw3D<T>::CirclePowerLaw3D(const SuperGeometry<T,3>& superGeometry,
188 int material, T maxVelocity, T n, T distance2Wall, T scale)
189 : AnalyticalF3D<T,T>(3), _maxVelocity(maxVelocity), _n(n), _scale(scale)
190{
191 _center = superGeometry.getStatistics().getCenterPhysR(material);
192 std::vector<T> normalTmp;
193 normalTmp.push_back(superGeometry.getStatistics().computeDiscreteNormal(material)[0]);
194 normalTmp.push_back(superGeometry.getStatistics().computeDiscreteNormal(material)[1]);
195 normalTmp.push_back(superGeometry.getStatistics().computeDiscreteNormal(material)[2]);
196 _normal = util::normalize(normalTmp);
197
198 // div. by 2 since one of the discrete redius directions is always 0 while the two other are not
199 _radius = T(distance2Wall);
200 for (int iD = 0; iD < 3; iD++) {
201 _radius += superGeometry.getStatistics().getPhysRadius(material)[iD]/2.;
202 }
203}
204
205template <typename T>
206CirclePowerLaw3D<T>::CirclePowerLaw3D(bool useMeanVelocity, std::vector<T> axisPoint, std::vector<T> axisDirection, T Velocity, T radius, T n, T scale)
207 : CirclePowerLaw3D<T>(axisPoint, axisDirection, Velocity, radius, n, scale)
208{
209 if (useMeanVelocity) {
210 _maxVelocity = (2. + (1. + n)/n)/((1. + n)/n * util::pow(1.,(2. + (1. + n)/n))) * Velocity;
211 }
212}
213
214template <typename T>
215CirclePowerLaw3D<T>::CirclePowerLaw3D(bool useMeanVelocity, T center0, T center1, T center2, T normal0, T normal1, T normal2, T radius, T Velocity, T n, T scale)
216 : CirclePowerLaw3D<T>(center0, center1, center2, normal0, normal1, normal2, radius, Velocity, n, scale)
217{
218 if (useMeanVelocity) {
219 _maxVelocity = (2. + (1. + n)/n)/((1. + n)/n * util::pow(1.,(2. + (1. + n)/n))) * Velocity;
220 }
221}
222
223template <typename T>
224CirclePowerLaw3D<T>::CirclePowerLaw3D(bool useMeanVelocity, SuperGeometry<T,3>& superGeometry, int material, T Velocity, T n, T distance2Wall, T scale)
225 : CirclePowerLaw3D<T>(superGeometry, material, Velocity, n, distance2Wall, scale)
226{
227 if (useMeanVelocity) {
228 _maxVelocity = (2. + (1. + n)/n)/((1. + n)/n * util::pow(1.,(2. + (1. + n)/n))) * Velocity;
229 }
230}
231
232template <typename T>
233bool CirclePowerLaw3D<T>::operator()(T output[], const T x[])
234{
235 output[0] = _scale*_maxVelocity*_normal[0]*(1.-util::pow(util::sqrt((x[1]-_center[1])*(x[1]-_center[1])+(x[2]-_center[2])*(x[2]-_center[2]))/_radius, (_n + 1.)/_n));
236 output[1] = _scale*_maxVelocity*_normal[1]*(1.-util::pow(util::sqrt((x[0]-_center[0])*(x[0]-_center[0])+(x[2]-_center[2])*(x[2]-_center[2]))/_radius, (_n + 1.)/_n));
237 output[2] = _scale*_maxVelocity*_normal[2]*(1.-util::pow(util::sqrt((x[1]-_center[1])*(x[1]-_center[1])+(x[0]-_center[0])*(x[0]-_center[0]))/_radius, (_n + 1.)/_n));
238 return true;
239}
240
241template <typename T>
242CirclePowerLawTurbulent3D<T>::CirclePowerLawTurbulent3D(std::vector<T> center, std::vector<T> normal,
243 T maxVelocity, T radius, T n, T turbulenceIntensity, T scale)
244 : CirclePowerLaw3D<T>(center, normal, maxVelocity, radius, n, scale), _turbulenceIntensity(turbulenceIntensity), _generator(_rd()), _dist(0.0, 1.0) { }
245
246template <typename T>
247CirclePowerLawTurbulent3D<T>::CirclePowerLawTurbulent3D(T center0, T center1, T center2, T normal0,
248 T normal1, T normal2, T radius, T maxVelocity, T n, T turbulenceIntensity, T scale)
249 : CirclePowerLaw3D<T>(center0, center1, center2, normal0, normal1, normal2, radius, maxVelocity, n, scale), _turbulenceIntensity(turbulenceIntensity), _generator(_rd()), _dist(0.0, 1.0) { }
250
251template <typename T>
252CirclePowerLawTurbulent3D<T>::CirclePowerLawTurbulent3D(SuperGeometry<T,3>& superGeometry,
253 int material, T maxVelocity, T distance2Wall, T n, T turbulenceIntensity, T scale)
254 : CirclePowerLaw3D<T>(superGeometry, material, maxVelocity, n, distance2Wall, scale), _turbulenceIntensity(turbulenceIntensity), _generator(_rd()), _dist(0.0, 1.0) { }
255
256template <typename T>
257CirclePowerLawTurbulent3D<T>::CirclePowerLawTurbulent3D(bool useMeanVelocity, std::vector<T> axisPoint, std::vector<T> axisDirection, T Velocity, T radius, T n, T turbulenceIntensity, T scale)
258 : CirclePowerLawTurbulent3D(axisPoint, axisDirection, Velocity, radius, n, turbulenceIntensity, scale)
259{
260 if (useMeanVelocity) {
261 this->_maxVelocity = ((1. + 1./n) * (2. + 1./n)) / 2. * Velocity;
262 }
263}
264template <typename T>
265CirclePowerLawTurbulent3D<T>::CirclePowerLawTurbulent3D(bool useMeanVelocity, T center0, T center1, T center2, T normal0, T normal1, T normal2, T radius, T Velocity, T n, T turbulenceIntensity, T scale)
266 : CirclePowerLawTurbulent3D(center0, center1, center2, normal0, normal1, normal2, radius, Velocity, n, turbulenceIntensity, scale)
267{
268 if (useMeanVelocity) {
269 this->_maxVelocity = ((1. + 1./n) * (2. + 1./n)) / 2. * Velocity;
270 }
271}
272
273template <typename T>
274CirclePowerLawTurbulent3D<T>::CirclePowerLawTurbulent3D(bool useMeanVelocity, SuperGeometry<T,3>& superGeometry, int material, T Velocity, T distance2Wall, T n, T turbulenceIntensity, T scale)
275 : CirclePowerLawTurbulent3D(superGeometry, material, Velocity, n, turbulenceIntensity, distance2Wall, scale)
276{
277 if (useMeanVelocity) {
278 this->_maxVelocity = ((1. + 1./n) * (2. + 1./n)) / 2. * Velocity;
279 }
280}
281
282template <typename T>
283bool CirclePowerLawTurbulent3D<T>::operator()(T output[], const T x[])
284{
285 T meanVelocity = this->_maxVelocity * 2.0 / ((1. + 1./this->_n) * (2. + 1./this->_n));
286
287 if ( 1.-util::sqrt((x[1]-this->_center[1])*(x[1]-this->_center[1])+(x[2]-this->_center[2])*(x[2]-this->_center[2]))/this->_radius < 0.) {
288 output[0] = T();
289 }
290 else {
291 output[0] = this->_scale*this->_maxVelocity*this->_normal[0]*
292 (util::pow(1.-util::sqrt((x[1]-this->_center[1])*(x[1]-this->_center[1])+(x[2]-this->_center[2])*(x[2]-this->_center[2]))/this->_radius, 1./this->_n));
293 output[0] += _dist(_generator) * _turbulenceIntensity * meanVelocity;
294 }
295 if ( 1.-util::sqrt((x[0]-this->_center[0])*(x[0]-this->_center[0])+(x[2]-this->_center[2])*(x[2]-this->_center[2]))/this->_radius < 0.) {
296 output[1] = T();
297 }
298 else {
299 output[1] = this->_scale*this->_maxVelocity*this->_normal[1]*
300 (util::pow(1.-util::sqrt((x[0]-this->_center[0])*(x[0]-this->_center[0])+(x[2]-this->_center[2])*(x[2]-this->_center[2]))/this->_radius, 1./this->_n));
301 output[1] += _dist(_generator) * _turbulenceIntensity * meanVelocity;
302 }
303 if ( 1.-util::sqrt((x[1]-this->_center[1])*(x[1]-this->_center[1])+(x[0]-this->_center[0])*(x[0]-this->_center[0]))/this->_radius < 0.) {
304 output[2] = T();
305 }
306 else {
307 output[2] = this->_scale*this->_maxVelocity*this->_normal[2]*
308 (util::pow(1.-util::sqrt((x[1]-this->_center[1])*(x[1]-this->_center[1])+(x[0]-this->_center[0])*(x[0]-this->_center[0]))/this->_radius, 1./this->_n));
309 output[2] += _dist(_generator) * _turbulenceIntensity * meanVelocity;
310 }
311
312 return true;
313}
314
315
316template <typename T>
317CirclePoiseuille3D<T>::CirclePoiseuille3D(std::vector<T> center, std::vector<T> normal,
318 T maxVelocity, T radius, T scale)
319 : CirclePowerLaw3D<T>(center, normal, maxVelocity, radius, 1., scale) { }
320
321template <typename T>
322CirclePoiseuille3D<T>::CirclePoiseuille3D(T center0, T center1, T center2, T normal0,
323 T normal1, T normal2, T radius, T maxVelocity, T scale)
324 : CirclePowerLaw3D<T>(center0, center1, center2, normal0, normal1, normal2, radius, maxVelocity, 1., scale) { }
325
326template <typename T>
327CirclePoiseuille3D<T>::CirclePoiseuille3D(SuperGeometry<T,3>& superGeometry,
328 int material, T maxVelocity, T distance2Wall, T scale)
329 : CirclePowerLaw3D<T>(superGeometry, material, maxVelocity, 1., distance2Wall, scale) { }
330
331template <typename T>
332CirclePoiseuille3D<T>::CirclePoiseuille3D(bool useMeanVelocity, std::vector<T> axisPoint, std::vector<T> axisDirection, T Velocity, T radius, T scale)
333 : CirclePoiseuille3D(axisPoint, axisDirection, Velocity, radius, scale)
334{
335 if (useMeanVelocity) {
336 this->_maxVelocity = 2. * Velocity;
337 }
338}
339
340template <typename T>
341CirclePoiseuille3D<T>::CirclePoiseuille3D(bool useMeanVelocity, T center0, T center1, T center2, T normal0, T normal1, T normal2, T radius, T Velocity, T scale)
342 : CirclePoiseuille3D(center0, center1, center2, normal0, normal1, normal2, radius, Velocity, scale)
343{
344 if (useMeanVelocity) {
345 this->_maxVelocity = 2. * Velocity;
346 }
347}
348
349template <typename T>
350CirclePoiseuille3D<T>::CirclePoiseuille3D(bool useMeanVelocity, SuperGeometry<T,3>& superGeometry, int material, T Velocity, T distance2Wall, T scale)
351 : CirclePoiseuille3D(superGeometry, material, Velocity, distance2Wall, scale)
352{
353 if (useMeanVelocity) {
354 this->_maxVelocity = 2. * Velocity;
355 }
356}
357
358template < typename T,typename DESCRIPTOR>
359CirclePoiseuilleStrainRate3D<T, DESCRIPTOR>::CirclePoiseuilleStrainRate3D(UnitConverter<T, DESCRIPTOR> const& converter, T ly) : AnalyticalF3D<T,T>(9)
360{
361 lengthY = ly;
362 lengthZ = ly;
363 maxVelocity = converter.getCharPhysVelocity();
364 this->getName() = "CirclePoiseuilleStrainRate3D";
365}
366
367
368template < typename T,typename DESCRIPTOR>
369bool CirclePoiseuilleStrainRate3D<T, DESCRIPTOR>::operator()(T output[], const T input[])
370{
371 T y = input[1];
372 T z = input[2];
373
374 T DuDx = T();
375 T DuDy = (T) maxVelocity*(-2.*(y-(lengthY))/((lengthY)*(lengthY)));
376 T DuDz = (T) maxVelocity*(-2.*(z-(lengthZ))/((lengthZ)*(lengthZ)));
377 T DvDx = T();
378 T DvDy = T();
379 T DvDz = T();
380 T DwDx = T();
381 T DwDy = T();
382 T DwDz = T();
383
384 output[0] = (DuDx + DuDx)/2.;
385 output[1] = (DuDy + DvDx)/2.;
386 output[2] = (DuDz + DwDx)/2.;
387 output[3] = (DvDx + DuDy)/2.;
388 output[4] = (DvDy + DvDy)/2.;
389 output[5] = (DvDz + DwDy)/2.;
390 output[6] = (DwDx + DuDz)/2.;
391 output[7] = (DwDy + DvDz)/2.;
392 output[8] = (DwDz + DwDz)/2.;
393
394 return true;
395};
396
397
398
399
400template <typename T>
401RectanglePoiseuille3D<T>::RectanglePoiseuille3D(olb::Vector<T, 3>& x0_, olb::Vector<T, 3>& x1_,
402 olb::Vector<T, 3>& x2_, std::vector<T>& maxVelocity_)
403 : AnalyticalF3D<T,T>(3), clout(std::cout,"RectanglePoiseille3D"), x0(x0_), x1(x1_),
404 x2(x2_), maxVelocity(maxVelocity_) { }
405
406
407template <typename T>
408RectanglePoiseuille3D<T>::RectanglePoiseuille3D(SuperGeometry<T,3>& superGeometry_,
409 int material_, std::vector<T>& maxVelocity_, T offsetX, T offsetY, T offsetZ)
410 : AnalyticalF3D<T,T>(3), clout(std::cout, "RectanglePoiseuille3D"), maxVelocity(maxVelocity_)
411{
412 olb::Vector<T, 3> min = superGeometry_.getStatistics().getMinPhysR(material_);
413 olb::Vector<T, 3> max = superGeometry_.getStatistics().getMaxPhysR(material_);
414 // set three corners of the plane, move by offset to land on the v=0
415 // boundary and adapt certain values depending on the orthogonal axis to get
416 // different points
417 x0 = min;
418 x1 = min;
419 x2 = min;
420 if ( util::nearZero(max[0]-min[0]) ) { // orthogonal to x-axis
421 x0[1] -= offsetY;
422 x0[2] -= offsetZ;
423 x1[1] -= offsetY;
424 x1[2] -= offsetZ;
425 x2[1] -= offsetY;
426 x2[2] -= offsetZ;
427
428 x1[1] = max[1] + offsetY;
429 x2[2] = max[2] + offsetZ;
430 }
431 else if ( util::nearZero(max[1]-min[1]) ) { // orthogonal to y-axis
432 x0[0] -= offsetX;
433 x0[2] -= offsetZ;
434 x1[0] -= offsetX;
435 x1[2] -= offsetZ;
436 x2[0] -= offsetX;
437 x2[2] -= offsetZ;
438
439 x1[0] = max[0] + offsetX;
440 x2[2] = max[2] + offsetZ;
441 }
442 else if ( util::nearZero(max[2]-min[2]) ) { // orthogonal to z-axis
443 x0[0] -= offsetX;
444 x0[1] -= offsetY;
445 x1[0] -= offsetX;
446 x1[1] -= offsetY;
447 x2[0] -= offsetX;
448 x2[1] -= offsetY;
449
450 x1[0] = max[0] + offsetX;
451 x2[1] = max[1] + offsetY;
452 }
453 else {
454 clout << "Error: constructor from material number works only for axis-orthogonal planes" << std::endl;
455 }
456}
457
458
459template <typename T>
460bool RectanglePoiseuille3D<T>::operator()(T output[], const T x[])
461{
462 std::vector<T> velocity(3,T());
463
464 // create plane normals and supports, (E1,E2) and (E3,E4) are parallel
465 std::vector<std::vector<T> > n(4,std::vector<T>(3,0)); // normal vectors
466 std::vector<std::vector<T> > s(4,std::vector<T>(3,0)); // support vectors
467 for (int iD = 0; iD < 3; iD++) {
468 n[0][iD] = x1[iD] - x0[iD];
469 n[1][iD] = -x1[iD] + x0[iD];
470 n[2][iD] = x2[iD] - x0[iD];
471 n[3][iD] = -x2[iD] + x0[iD];
472 s[0][iD] = x0[iD];
473 s[1][iD] = x1[iD];
474 s[2][iD] = x0[iD];
475 s[3][iD] = x2[iD];
476 }
477 for (int iP = 0; iP < 4; iP++) {
478 n[iP] = util::normalize(n[iP]);
479 }
480
481 // determine plane coefficients in the coordinate equation E_i: Ax+By+Cz+D=0
482 // given form: (x-s)*n=0 => A=n1, B=n2, C=n3, D=-(s1n1+s2n2+s3n3)
483 std::vector<T> A(4,0);
484 std::vector<T> B(4,0);
485 std::vector<T> C(4,0);
486 std::vector<T> D(4,0);
487
488 for (int iP = 0; iP < 4; iP++) {
489 A[iP] = n[iP][0];
490 B[iP] = n[iP][1];
491 C[iP] = n[iP][2];
492 D[iP] = -(s[iP][0]*n[iP][0] + s[iP][1]*n[iP][1] + s[iP][2]*n[iP][2]);
493 }
494
495 // determine distance to the walls by formula
496 std::vector<T> d(4,0);
497 T aabbcc(0);
498 for (int iP = 0; iP < 4; iP++) {
499 aabbcc = A[iP]*A[iP] + B[iP]*B[iP] + C[iP]*C[iP];
500 d[iP] = util::fabs(A[iP]*x[0]+B[iP]*x[1]+C[iP]*x[2]+D[iP])*util::pow(aabbcc,-0.5);
501 }
502
503 // length and width of the rectangle
504 std::vector<T> l(2,0);
505 l[0] = d[0] + d[1];
506 l[1] = d[2] + d[3];
507
508 T positionFactor = 16.*d[0]/l[0]*d[1]/l[0]*d[2]/l[1]*d[3]/l[1]; // between 0 and 1
509
510 output[0] = maxVelocity[0]*positionFactor;
511 output[1] = maxVelocity[1]*positionFactor;
512 output[2] = maxVelocity[2]*positionFactor;
513 return true;
514}
515
516//Trigonometric version TODO: to be unified( etc. plane calculation )
517template <typename T>
518RectangleTrigonometricPoiseuille3D<T>::RectangleTrigonometricPoiseuille3D(std::vector<T>& x0_, std::vector<T>& x1_,
519 std::vector<T>& x2_, std::vector<T>& maxVelocity_, int numOfPolynoms_)
520 : AnalyticalF3D<T,T>(3), clout(std::cout,"RectangleTrigonometricPoiseille3D"), x0(x0_), x1(x1_),
521 x2(x2_), maxVelocity(maxVelocity_), numOfPolynoms(numOfPolynoms_)
522{
523 calcPeakAndAvg();
524}
525
526
527template <typename T>
528RectangleTrigonometricPoiseuille3D<T>::RectangleTrigonometricPoiseuille3D(SuperGeometry<T,3>& superGeometry_,
529 int material_, std::vector<T>& maxVelocity_, T offsetX, T offsetY, T offsetZ, int numOfPolynoms_)
530 : AnalyticalF3D<T,T>(3), clout(std::cout, "RectangleTrigonometricPoiseuille3D"), maxVelocity(maxVelocity_),
531 numOfPolynoms(numOfPolynoms_)
532{
533 std::vector<T> min = superGeometry_.getStatistics().getMinPhysR(material_);
534 std::vector<T> max = superGeometry_.getStatistics().getMaxPhysR(material_);
535 // set three corners of the plane, move by offset to land on the v=0
536 // boundary and adapt certain values depending on the orthogonal axis to get
537 // different points
538 x0 = min;
539 x1 = min;
540 x2 = min;
541 if ( util::nearZero(max[0]-min[0]) ) { // orthogonal to x-axis
542 x0[1] -= offsetY;
543 x0[2] -= offsetZ;
544 x1[1] -= offsetY;
545 x1[2] -= offsetZ;
546 x2[1] -= offsetY;
547 x2[2] -= offsetZ;
548
549 x1[1] = max[1] + offsetY;
550 x2[2] = max[2] + offsetZ;
551 }
552 else if ( util::nearZero(max[1]-min[1]) ) { // orthogonal to y-axis
553 x0[0] -= offsetX;
554 x0[2] -= offsetZ;
555 x1[0] -= offsetX;
556 x1[2] -= offsetZ;
557 x2[0] -= offsetX;
558 x2[2] -= offsetZ;
559
560 x1[0] = max[0] + offsetX;
561 x2[2] = max[2] + offsetZ;
562 }
563 else if ( util::nearZero(max[2]-min[2]) ) { // orthogonal to z-axis
564 x0[0] -= offsetX;
565 x0[1] -= offsetY;
566 x1[0] -= offsetX;
567 x1[1] -= offsetY;
568 x2[0] -= offsetX;
569 x2[1] -= offsetY;
570
571 x1[0] = max[0] + offsetX;
572 x2[1] = max[1] + offsetY;
573 }
574 else {
575 clout << "Error: constructor from material number works only for axis-orthogonal planes" << std::endl;
576 }
577 calcPeakAndAvg();
578}
579
580
581template <typename T>
582void RectangleTrigonometricPoiseuille3D<T>::calcPeakAndAvg()
583{
584 //2006_Bruus_Lecture Theoretical microfluidics
585 T sumMax = 0;
586 T sumAvg = 0;
587 for (int i=1; i<=numOfPolynoms; ++i) { //TODO: higher polynoms still to be tested
588 T n = i*2-1; //Creating only odd indices
589 T A = util::sin(n*M_PI*0.5);
590 T max = (1./util::pow(n,3)) * (1.-(1 / util::cosh(n*M_PI / 2 ))) * A;
591 //TODO: only for y,z=1 //TODO: Maximum calculation needs to be adapted for different edge ratios
592 T avg = 2. * A * A * (n*M_PI-2*util::tanh(n*M_PI*0.5)) / (util::pow(n,5)*M_PI*M_PI);
593 sumMax += max;
594 sumAvg += avg;
595 }
596 profilePeak = sumMax;
597 profileRatioAvgToPeak = sumAvg / sumMax;
598}
599
600
601template <typename T>
602T RectangleTrigonometricPoiseuille3D<T>::getProfilePeak()
603{
604 return profilePeak;
605}
606
607template <typename T>
608T RectangleTrigonometricPoiseuille3D<T>::getRatioAvgToPeak()
609{
610 return profileRatioAvgToPeak;
611}
612
613template <typename T> //TODO: for now only works in x direction!
614bool RectangleTrigonometricPoiseuille3D<T>::operator()(T output[], const T x[])
615{
616 OstreamManager clout( std::cout,"RectPois" ); //TODO: to be removed
617
618 // create plane normals and supports, (E1,E2) and (E3,E4) are parallel
619 std::vector<std::vector<T> > n(4,std::vector<T>(3,0)); // normal vectors
620 std::vector<std::vector<T> > s(4,std::vector<T>(3,0)); // support vectors
621 for (int iD = 0; iD < 3; iD++) {
622 n[0][iD] = x1[iD] - x0[iD];
623 n[1][iD] = -x1[iD] + x0[iD];
624 n[2][iD] = x2[iD] - x0[iD];
625 n[3][iD] = -x2[iD] + x0[iD];
626 s[0][iD] = x0[iD];
627 s[1][iD] = x1[iD];
628 s[2][iD] = x0[iD];
629 s[3][iD] = x2[iD];
630 }
631 for (int iP = 0; iP < 4; iP++) {
632 n[iP] = util::normalize(n[iP]);
633 }
634
635 // determine plane coefficients in the coordinate equation E_i: Ax+By+Cz+D=0
636 // given form: (x-s)*n=0 => A=n1, B=n2, C=n3, D=-(s1n1+s2n2+s3n3)
637 std::vector<T> A(4,0);
638 std::vector<T> B(4,0);
639 std::vector<T> C(4,0);
640 std::vector<T> D(4,0);
641
642 for (int iP = 0; iP < 4; iP++) {
643 A[iP] = n[iP][0];
644 B[iP] = n[iP][1];
645 C[iP] = n[iP][2];
646 D[iP] = -(s[iP][0]*n[iP][0] + s[iP][1]*n[iP][1] + s[iP][2]*n[iP][2]);
647 }
648
649 // determine distance to the walls by formula
650 std::vector<T> d(4,0);
651 T aabbcc(0);
652 for (int iP = 0; iP < 4; iP++) {
653 aabbcc = A[iP]*A[iP] + B[iP]*B[iP] + C[iP]*C[iP];
654 d[iP] = util::fabs(A[iP]*x[0]+B[iP]*x[1]+C[iP]*x[2]+D[iP])*util::pow(aabbcc,-0.5);
655 }
656
657 // length and width of the rectangle
658 std::vector<T> l(2,0);
659 l[0] = d[0] + d[1] ;
660 l[1] = d[2] + d[3] ;
661
662 std::vector<T> dl(2,0);
663 dl[0] = d[0] / l[0];
664 dl[1] = d[2] / l[1];
665
666 //Polynomial creation for odd indices (2i-i) from:
667 //2006_Bruus_Lecture Theoretical microfluidics
668 T sumPoly = 0;
669 for (int i=1; i<=numOfPolynoms; ++i) { //TODO: higher polynoms still to be tested
670 T n = i*2-1;
671 T poly = (1./util::pow(n,3)) * (1.-(util::cosh(n*M_PI*(dl[0]-0.5)) / util::cosh(n*M_PI/ (2) ))) * util::sin(n*M_PI*dl[1]);
672 sumPoly += poly;
673 }
674
675 output[0] = maxVelocity[0] * ( sumPoly / profilePeak);
676 output[1] = maxVelocity[1] * ( sumPoly / profilePeak);
677 output[2] = maxVelocity[2] * ( sumPoly / profilePeak);
678 return true;
679}
680
681
682template <typename T>
683EllipticPoiseuille3D<T>::EllipticPoiseuille3D(std::vector<T> center, T a, T b, T maxVel)
684 : AnalyticalF3D<T,T>(3), clout(std::cout, "EllipticPoiseuille3D"), _center(center),
685 _a2(a*a), _b2(b*b), _maxVel(maxVel) { }
686
687
688template <typename T>
689bool EllipticPoiseuille3D<T>::operator()(T output[], const T x[])
690{
691 T cosOmega2 = util::pow(x[0]-_center[0],2.)/(util::pow(x[0]-_center[0],2.)+util::pow(x[1]-_center[1],2.));
692 T rad2 = _a2*_b2 /((_b2-_a2)*cosOmega2 + _a2) ;
693 T x2 = (util::pow(x[0]-_center[0],2.)+util::pow(x[1]-_center[1],2.))/rad2;
694
695 // clout << _center[0] << " " << _center[1] << " " << x[0] << " " << x[1] << " " << util::sqrt(rad2) << " " << util::sqrt(util::pow(x[0]-_center[0],2.)+util::pow(x[1]-_center[1],2.)) << " " << x2 << std::endl;
696
697 output[0] = 0.;
698 output[1] = 0.;
699 output[2] = _maxVel*(-x2+1);
700 if ( util::nearZero(_center[0]-x[0]) && util::nearZero(_center[1]-x[1]) ) {
701 output[2] = _maxVel;
702 }
703 return true;
704}
705
706
707template <typename T>
708AnalyticalPorousVelocity3D<T>::AnalyticalPorousVelocity3D(SuperGeometry<T,3>& superGeometry, int material, T K_, T mu_, T gradP_, T radius_, T eps_)
709 : AnalyticalF3D<T,T>(3), K(K_), mu(mu_), gradP(gradP_), radius(radius_), eps(eps_)
710{
711 this->getName() = "AnalyticalPorousVelocity3D";
712
713 center = superGeometry.getStatistics().getCenterPhysR(material);
714 std::vector<T> normalTmp;
715 normalTmp.push_back(superGeometry.getStatistics().computeDiscreteNormal(material)[0]);
716 normalTmp.push_back(superGeometry.getStatistics().computeDiscreteNormal(material)[1]);
717 normalTmp.push_back(superGeometry.getStatistics().computeDiscreteNormal(material)[2]);
718 normal = util::normalize(normalTmp);
719
720};
721
722
723template <typename T>
724T AnalyticalPorousVelocity3D<T>::getPeakVelocity()
725{
726 T uMax = K / mu*gradP*(1. - 1./(util::cosh((util::sqrt(1./K))*radius)));
727
728 return uMax/eps;
729};
730
731
732template <typename T>
733bool AnalyticalPorousVelocity3D<T>::operator()(T output[], const T x[])
734{
735 T dist[3] = {};
736 dist[0] = normal[0]*util::sqrt((x[1]-center[1])*(x[1]-center[1])+(x[2]-center[2])*(x[2]-center[2]));
737 dist[1] = normal[1]*util::sqrt((x[0]-center[0])*(x[0]-center[0])+(x[2]-center[2])*(x[2]-center[2]));
738 dist[2] = normal[2]*util::sqrt((x[1]-center[1])*(x[1]-center[1])+(x[0]-center[0])*(x[0]-center[0]));
739
740 output[0] = K / mu*gradP*(1. - (util::cosh((util::sqrt(1./K))*(dist[0])))/(util::cosh((util::sqrt(1./K))*radius)));
741 output[1] = K / mu*gradP*(1. - (util::cosh((util::sqrt(1./K))*(dist[1])))/(util::cosh((util::sqrt(1./K))*radius)));
742 output[2] = K / mu*gradP*(1. - (util::cosh((util::sqrt(1./K))*(dist[2])))/(util::cosh((util::sqrt(1./K))*radius)));
743
744 output[0] *= normal[0]/eps;
745 output[1] *= normal[1]/eps;
746 output[2] *= normal[2]/eps;
747
748 return true;
749};
750
751
752
754
755
756// constructor defines _referenceVector to obtain angle between 2 vectors,
757// vector x has to be turned by angle in mathematical positive sense depending
758// to _orientation to lie over _referenceVector
759template<typename T, typename S>
760AngleBetweenVectors3D<T, S>::AngleBetweenVectors3D(
761 std::vector<T> referenceVector, std::vector<T> orientation)
762 : AnalyticalF3D<T, S>(1),
763 _referenceVector(referenceVector),
764 _orientation(orientation)
765{
766 this->getName() = "const";
767}
768
769// operator returns angle between vectors x and _referenceVector
770template<typename T, typename S>
771bool AngleBetweenVectors3D<T, S>::operator()(T output[], const S x[])
772{
773 Vector<T, 3> n_x;
774 Vector<T, 3> orientation(_orientation);
775 T angle = T(0);
776
777 if ( util::nearZero(x[0]) && util::nearZero(x[1]) && util::nearZero(x[2]) ) {
778 // if (util::abs(x[0]) + util::abs(x[1]) + util::abs(x[2]) == T()) {
779 output[0] = angle; // angle = 0
780 return true;
781 }
782 else {
783 //Vector<S, 3> x_tmp(x); // check construction
784 n_x[0] = x[0];
785 n_x[1] = x[1];
786 n_x[2] = x[2];
787 n_x = normalize(n_x);
788 }
789
790 Vector<T, 3> n_ref(_referenceVector);
791 n_ref = normalize(n_ref);
792 Vector<T, 3> cross = crossProduct3D(n_x, n_ref);
793 T n_dot = n_x * n_ref;
794
795
796 if ( util::nearZero(cross*orientation) ) {
797 // std::cout<< "orientation in same plane as x and refvector" << std::endl;
798 }
799 // angle = Pi, if n_x, n_ref opposite
800 if ( util::nearZero(n_x[0]+n_ref[0]) && util::nearZero(n_x[1]+n_ref[1]) && util::nearZero(n_x[2]+n_ref[2]) ) {
801 angle = util::acos(-1);
802 }
803 // angle = 0, if n_x, n_ref equal
804 else if ( util::nearZero(n_x[0]-n_ref[0]) && util::nearZero(n_x[1]-n_ref[1]) && util::nearZero(n_x[2]-n_ref[2]) ) {
805 angle = T();
806 }
807 // angle in (0,Pi) or (Pi,2Pi), if n_x, n_ref not opposite or equal
808 else {
809 Vector<T, 3> n_cross(cross);
810 n_cross = normalize(n_cross);
811 T normal = norm(cross);
812 Vector<T, 3> n_orient;
813
814 if ( !util::nearZero(norm(orientation)) ) {
815 n_orient = orientation;
816 n_orient = normalize(n_orient);
817 }
818 else {
819 std::cout << "orientation vector does not fit" << std::endl;
820 }
821 if ((cross * orientation) > T()) {
822 angle = 2*M_PI - util::atan2(normal, n_dot);
823 }
824 if ((cross * orientation) < T()) {
825 angle = util::atan2(normal, n_dot);
826 }
827 }
828 output[0] = angle;
829 return true;
830}
831
832
833// constructor to obtain rotation util::round an _rotAxisDirection with angle _alpha
834// and _origin
835template<typename T, typename S>
836RotationRoundAxis3D<T, S>::RotationRoundAxis3D(olb::Vector<T, 3> origin,
837 olb::Vector<T, 3> rotAxisDirection, T alpha)
838 : AnalyticalF3D<T, S>(3),
839 _origin(origin),
840 _rotAxisDirection(rotAxisDirection),
841 _alpha(alpha)
842{
843 this->getName() = "const";
844}
845
846// operator returns coordinates of the resulting point after rotation of x
847// with _alpha in math positive sense to _rotAxisDirection through _origin
848template<typename T, typename S>
849bool RotationRoundAxis3D<T, S>::operator()(T output[], const S x[])
850{
851 Vector<T, 3> n(_rotAxisDirection);
852 if ( !util::nearZero(norm(n)) && norm(n) > 0 ) {
853 //std::cout<< "Rotation axis: " << _rotAxisDirection[0] << " " << _rotAxisDirection[1] << " " << _rotAxisDirection[2] << std::endl;
854 n = normalize(n);
855 // translation
856 Vector<T, 3> x_tmp;
857 for (int i = 0; i < 3; ++i) {
858 x_tmp[i] = x[i] - _origin[i];
859 }
860 // rotation
861 output[0] = (n[0]*n[0]*(1 - util::cos(_alpha)) + util::cos(_alpha)) * x_tmp[0]
862 + (n[0]*n[1]*(1 - util::cos(_alpha)) - n[2]*util::sin(_alpha))*x_tmp[1]
863 + (n[0]*n[2]*(1 - util::cos(_alpha)) + n[1]*util::sin(_alpha))*x_tmp[2]
864 + _origin[0];
865
866 output[1] = (n[1]*n[0]*(1 - util::cos(_alpha)) + n[2]*util::sin(_alpha))*x_tmp[0]
867 + (n[1]*n[1]*(1 - util::cos(_alpha)) + util::cos(_alpha))*x_tmp[1]
868 + (n[1]*n[2]*(1 - util::cos(_alpha)) - n[0]*util::sin(_alpha))*x_tmp[2]
869 + _origin[1];
870
871 output[2] = (n[2]*n[0]*(1 - util::cos(_alpha)) - n[1]*util::sin(_alpha))*x_tmp[0]
872 + (n[2]*n[1]*(1 - util::cos(_alpha)) + n[0]*util::sin(_alpha))*x_tmp[1]
873 + (n[2]*n[2]*(1 - util::cos(_alpha)) + util::cos(_alpha))*x_tmp[2]
874 + _origin[2];
875 return true;
876 }
877 else {
878 //std::cout<< "Axis is null or NaN" <<std::endl;
879 for (int j=0; j<3; j++) {
880 output[j] = x[j];
881 }
882 return true;
883 }
884}
885
886
888
889
890// constructor to obtain Cartesian coordinates of cylindrical coordinates,
891// the z-axis direction equals the axis direction of the cylinder
892// with _cylinderOrigin
893template<typename T, typename S>
894CylinderToCartesian3D<T, S>::CylinderToCartesian3D(
895 std::vector<T> cylinderOrigin)
896 : AnalyticalF3D<T, S>(3),
897 _cylinderOrigin(cylinderOrigin)
898{
899 this->getName() = "CylinderToCartesian3D";
900}
901
902// operator returns Cartesian coordinates of cylindrical coordinates,
903// input is x[0] = radius, x[1] = phi, x[2] = z
904template<typename T, typename S>
905bool CylinderToCartesian3D<T, S>::operator()(T output[], const S x[])
906{
907 PolarToCartesian2D<T, S> pol2cart(_cylinderOrigin);
908 pol2cart(output,x);
909 output[2] = x[2];
910 return true;
911}
912
913
914// constructor to obtain cylindrical coordinates of Cartesian coordinates
915template<typename T, typename S>
916CartesianToCylinder3D<T, S>::CartesianToCylinder3D(
917 olb::Vector<T, 3> cartesianOrigin, T& angle, olb::Vector<T, 3> orientation)
918 : AnalyticalF3D<T, S>(3),
919 _cartesianOrigin(cartesianOrigin),
920 _orientation(orientation)
921{
922 // rotation with angle to (0,0,1)
923 std::vector<T> origin(3,T());
924 T e3[3]= {T(),T(),T()};
925 e3[2] = T(1);
926 RotationRoundAxis3D<T, S> rotRAxis(origin, orientation, angle);
927 T tmp[3] = {T(),T(),T()};
928 rotRAxis(tmp, e3);
929
930 std::vector<T> htmp(tmp,tmp+3);
931 _axisDirection = htmp;
932}
933
934// constructor to obtain cylindrical coordinates of Cartesian coordinates
935template<typename T, typename S>
936CartesianToCylinder3D<T, S>::CartesianToCylinder3D(
937 olb::Vector<T, 3> cartesianOrigin, olb::Vector<T, 3> axisDirection,
938 olb::Vector<T, 3> orientation)
939 : AnalyticalF3D<T, S>(3),
940 _cartesianOrigin(cartesianOrigin),
941 _axisDirection(axisDirection),
942 _orientation(orientation)
943{
944 this->getName() = "const";
945}
946
947// constructor to obtain cylindrical coordinates of Cartesian coordinates
948template<typename T, typename S>
949CartesianToCylinder3D<T, S>::CartesianToCylinder3D(
950 T cartesianOriginX, T cartesianOriginY, T cartesianOriginZ,
951 T axisDirectionX, T axisDirectionY, T axisDirectionZ,
952 T orientationX, T orientationY, T orientationZ)
953 : AnalyticalF3D<T, S>(3)
954{
955 _cartesianOrigin[0] = cartesianOriginX;
956 _cartesianOrigin[1] = cartesianOriginY;
957 _cartesianOrigin[2] = cartesianOriginZ;
958
959 _axisDirection[0] = axisDirectionX;
960 _axisDirection[1] = axisDirectionY;
961 _axisDirection[2] = axisDirectionZ;
962
963 _orientation[0] = orientationX;
964 _orientation[1] = orientationY;
965 _orientation[2] = orientationZ;
966}
967
968// operator returns cylindrical coordinates, output[0] = radius(>= 0),
969// output[1] = phi in [0, 2Pi), output[2] = z
970template<typename T, typename S>
971bool CartesianToCylinder3D<T, S>::operator()(T output[], const S x[])
972{
973 T x_rot[3] = {T(),T(),T()};
974 x_rot[0] = x[0];
975 x_rot[1] = x[1];
976 x_rot[2] = x[2];
977// T x_rot[3] = {x[0], x[1], x[2]};
978 Vector<T, 3> e3(T(),T(), T(1));
979 Vector<T, 3> axisDirection(_axisDirection);
980 Vector<T, 3> orientation(_orientation);
981
982 Vector<T, 3> normal = crossProduct3D(axisDirection, e3);
983 Vector<T, 3> normalAxisDir(axisDirection);
984 normalAxisDir = normalize(normalAxisDir);
985
986 // if axis has to be rotated
987 if (!( util::nearZero(normalAxisDir[0]) && util::nearZero(normalAxisDir[1]) && util::nearZero(normalAxisDir[2]-1) )) {
988
989 if ( !util::nearZero(norm(orientation)) ) {
990 normal = orientation; // if _orientation = 0,0,0 -> segm. fault
991 }
992 AngleBetweenVectors3D<T,S> angle(util::fromVector3(e3), util::fromVector3(normal));
993 T tmp[3] = {T(),T(),T()};
994 tmp[0] = axisDirection[0];
995 tmp[1] = axisDirection[1];
996 tmp[2] = axisDirection[2];
997 T alpha[1] = {T()};
998 angle(alpha, tmp);
999
1000 // rotation with angle alpha to (0,0,1)
1001 RotationRoundAxis3D<T, S> rotRAxis(_cartesianOrigin, util::fromVector3(normal), -alpha[0]);
1002 rotRAxis(x_rot, x);
1003 }
1004 CartesianToPolar2D<T, S> car2pol(_cartesianOrigin, util::fromVector3(e3), _orientation);
1005 T x_rot_help[2] = {T(),T()};
1006 x_rot_help[0] = x_rot[0];
1007 x_rot_help[1] = x_rot[1];
1008
1009 T output_help[2] = {T(),T()};
1010 car2pol(output_help, x_rot_help);
1011
1012 output[0] = output_help[0];
1013 output[1] = output_help[1];
1014 output[2] = x_rot[2] - _cartesianOrigin[2];
1015 return true; // output[0] = radius, output[1] = phi, output[2] = z
1016}
1017
1018// Read only access to _axisDirection
1019template<typename T, typename S>
1020olb::Vector<T, 3> const& CartesianToCylinder3D<T, S>::getAxisDirection() const
1021{
1022 return _axisDirection;
1023}
1024
1025// Read and write access to _axisDirection
1026template<typename T, typename S>
1027olb::Vector<T, 3>& CartesianToCylinder3D<T, S>::getAxisDirection()
1028{
1029 return _axisDirection;
1030}
1031
1032// Read only access to _catesianOrigin
1033template<typename T, typename S>
1034olb::Vector<T, 3> const& CartesianToCylinder3D<T, S>::getCartesianOrigin() const
1035{
1036 return _cartesianOrigin;
1037}
1038
1039// Read and write access to _catesianOrigin
1040template<typename T, typename S>
1041olb::Vector<T, 3>& CartesianToCylinder3D<T, S>::getCartesianOrigin()
1042{
1043 return _cartesianOrigin;
1044}
1045
1046// Read and write access to _axisDirection using angle
1047template<typename T, typename S>
1048void CartesianToCylinder3D<T, S>::setAngle(T angle)
1049{
1050 std::vector<T> Null(3,T());
1051 T e3[3] = {T(),T(),T()};
1052 e3[2] = T(1);
1053
1054 RotationRoundAxis3D<T, S> rotRAxis(Null, _orientation, angle);
1055 T tmp[3] = {T(),T(),T()};
1056 rotRAxis(tmp,e3);
1057 for (int i=0; i<3; ++i) {
1058 _axisDirection[i] = tmp[i];
1059 }
1060}
1061
1062
1064
1065
1066// constructor to obtain spherical coordinates of Cartesian coordinates
1067template<typename T, typename S>
1068SphericalToCartesian3D<T, S>::SphericalToCartesian3D()
1069//std::vector<T> sphericalOrigin)
1070 : AnalyticalF3D<T, S>(3) //, _sphericalOrigin(sphericalOrigin)
1071{
1072 this->getName() = "const";
1073}
1074
1075// operator returns Cartesian coordinates of spherical coordinates
1076// (x[0] = radius, x[1] = phi, x[2] = theta)
1077// phi in x-y-plane, theta in x-z-plane, z axis as orientation
1078template<typename T, typename S>
1079bool SphericalToCartesian3D<T, S>::operator()(T output[], const S x[])
1080{
1081 output[0] = x[0]*util::sin(x[2])*util::cos(x[1]);
1082 output[1] = x[0]*util::sin(x[2])*util::sin(x[1]);
1083 output[2] = x[0]*util::cos(x[2]);
1084 return true;
1085}
1086
1087
1088template<typename T, typename S>
1089CartesianToSpherical3D<T, S>::CartesianToSpherical3D(
1090 std::vector<T> cartesianOrigin, std::vector<T> axisDirection)
1091 : AnalyticalF3D<T, S>(3),
1092 _cartesianOrigin(cartesianOrigin),
1093 _axisDirection(axisDirection)
1094{
1095 this->getName() = "const";
1096}
1097
1098// operator returns spherical coordinates output[0] = radius(>= 0),
1099// output[1] = phi in [0, 2Pi), output[2] = theta in [0, Pi]
1100template<typename T, typename S>
1101bool CartesianToSpherical3D<T, S>::operator()(T output[], const S x[])
1102{
1103 T x_rot[3] = {T(),T(),T()};
1104 x_rot[0] = x[0];
1105 x_rot[1] = x[1];
1106 x_rot[2] = x[2];
1107
1108 Vector<T,3> axisDirection(_axisDirection);
1109 Vector<T,3> e3(T(), T(), T(1));
1110 Vector<T,3> normal = crossProduct3D(axisDirection,e3);
1111
1112 Vector<T,3> normalAxisDir = axisDirection;
1113 normalAxisDir = normalize(normalAxisDir);
1114 Vector<T,3> cross;
1115 // if axis has to be rotated
1116 if ( !( util::nearZero(normalAxisDir[0]) && util::nearZero(normalAxisDir[1]) && util::nearZero(normalAxisDir[2]-1) ) ) {
1117 AngleBetweenVectors3D<T, S> angle(util::fromVector3(e3), util::fromVector3(normal));
1118 T tmp[3] = {T(),T(),T()};
1119 for (int i=0; i<3; ++i) {
1120 tmp[i] = axisDirection[i];
1121 }
1122 T alpha[1] = {T(0)};
1123 angle(alpha,tmp);
1124 // cross is rotation axis
1125 cross = crossProduct3D(e3, axisDirection);
1126 // rotation with angle alpha to (0,0,1)
1127 RotationRoundAxis3D<T, S> rotRAxis(_cartesianOrigin, util::fromVector3(normal), -alpha[0]);
1128 rotRAxis(x_rot,x);
1129 }
1130
1131 CartesianToPolar2D<T, S> car2pol(_cartesianOrigin, util::fromVector3(e3), util::fromVector3(cross));
1132 // y[0] = car2pol(x_rot)[0];
1133 Vector<T, 3> difference;
1134
1135 for (int i=0; i<3; ++i) {
1136 difference[i] = x_rot[i] - _cartesianOrigin[i];
1137 }
1138
1139 T distance = T();
1140 for (int i = 0; i <= 2; ++i) {
1141 distance += difference[i]*difference[i];
1142 }
1143 distance = util::sqrt(distance);
1144
1145 car2pol(output, x_rot);
1146 output[0] = distance;
1147 output[2] = util::acos(difference[2]/output[0]); // angle between z-axis and r-vector
1148
1149 return true; // output[0] = radius, output[1] = phi, output[2] = theta
1150}
1151
1152
1154
1155
1161template<typename T, typename S>
1162MagneticForceFromCylinder3D<T, S>::MagneticForceFromCylinder3D(
1163 CartesianToCylinder3D<T, S>& car2cyl, T length, T radWire, T cutoff, T Mw, T Mp,
1164 T radParticle, T mu0)
1165 : AnalyticalF3D<T, S>(3),
1166 _car2cyl(car2cyl),
1167 _length(length),
1168 _radWire(radWire),
1169 _cutoff(cutoff),
1170 _Mw(Mw),
1171 _Mp(Mp),
1172 _radParticle(radParticle),
1173 _mu0(mu0)
1174{
1175 // Field direction H_0 parallel to fluid flow direction, if not: *(-1)
1176 _factor = -_mu0*4/3*M_PI*util::pow(_radParticle, 3)*_Mp*_Mw*util::pow(_radWire, 2);
1177}
1178
1179template<typename T, typename S>
1180bool MagneticForceFromCylinder3D<T, S>::operator()(T output[], const S x[])
1181{
1182
1183 Vector<T, 3> magneticForcePolar;
1184 T outputTmp[3] = {T(), T(), T()};
1185 Vector<T, 3> normalAxis(_car2cyl.getAxisDirection());
1186 normalAxis = normalize(normalAxis);
1187
1188 Vector<T, 3> relPosition;
1189 relPosition[0] = (x[0] - _car2cyl.getCartesianOrigin()[0]);
1190 relPosition[1] = (x[1] - _car2cyl.getCartesianOrigin()[1]);
1191 relPosition[2] = (x[2] - _car2cyl.getCartesianOrigin()[2]);
1192
1193 T tmp[3] = {T(),T(),T()};
1194 _car2cyl(tmp,x);
1195 T rad = tmp[0];
1196 T phi = tmp[1];
1197
1198 T test = relPosition * normalAxis;
1199
1200 if ( (rad > _radWire && rad <= _cutoff*_radWire) &&
1201 (T(0) <= test && test <= _length) ) {
1202
1203 magneticForcePolar[0] = _factor/util::pow(rad, 3)
1204 *(util::pow(_radWire, 2)/util::pow(rad, 2) + util::cos(2*phi));
1205 magneticForcePolar[1] = _factor/util::pow(rad, 3)*util::sin(2*phi);
1206
1207 // changes radial and angular force components to cartesian components.
1208 outputTmp[0] = magneticForcePolar[0]*util::cos(phi)
1209 - magneticForcePolar[1]*util::sin(phi);
1210 outputTmp[1] = magneticForcePolar[0]*util::sin(phi)
1211 + magneticForcePolar[1]*util::cos(phi);
1212
1213 // if not in standard axis direction
1214 if ( !( util::nearZero(normalAxis[0]) && util::nearZero(normalAxis[1]) && util::nearZero(normalAxis[2]-1) ) ) {
1215 Vector<T, 3> e3(T(), T(), T(1));
1216 Vector<T, 3> orientation =
1217 crossProduct3D(Vector<T,3>(_car2cyl.getAxisDirection()),e3);
1218
1219 AngleBetweenVectors3D<T,S> angle(util::fromVector3(e3), util::fromVector3(orientation));
1220 T alpha[1] = {T()};
1221 T tmp2[3] = {T(),T(),T()};
1222 for (int i = 0; i<3; ++i) {
1223 tmp2[i] = _car2cyl.getAxisDirection()[i];
1224 }
1225 angle(alpha,tmp2);
1226
1227 std::vector<T> origin(3, T());
1228 RotationRoundAxis3D<T, S> rotRAxis(origin, util::fromVector3(orientation), alpha[0]);
1229 rotRAxis(output,outputTmp);
1230 }
1231 else {
1232 output[0] = outputTmp[0];
1233 output[1] = outputTmp[1];
1234 output[2] = T();
1235 }
1236 }
1237 else {
1238 output[0] = outputTmp[0];
1239 output[1] = outputTmp[1];
1240 output[2] = outputTmp[1];
1241 }
1242 return true;
1243}
1244
1245template<typename T, typename S>
1246MagneticFieldFromCylinder3D<T, S>::MagneticFieldFromCylinder3D(
1247 CartesianToCylinder3D<T, S>& car2cyl, T length, T radWire, T cutoff, T Mw
1248)
1249 : AnalyticalF3D<T, S>(3),
1250 _car2cyl(car2cyl),
1251 _length(length),
1252 _radWire(radWire),
1253 _cutoff(cutoff),
1254 _Mw(Mw)
1255{
1256 // Field direction H_0 parallel to fluid flow direction, if not: *(-1)
1257 _factor = - _Mw * util::pow(_radWire, 2);
1258}
1259
1260template<typename T, typename S>
1261bool MagneticFieldFromCylinder3D<T, S>::operator()(T output[], const S x[])
1262{
1263
1264 Vector<T, 3> magneticFieldPolar;
1265 T outputTmp[3] = {T(), T(), T()};
1266 Vector<T, 3> normalAxis(_car2cyl.getAxisDirection());
1267 normalAxis = normalize(normalAxis);
1268
1269 Vector<T, 3> relPosition;
1270 relPosition[0] = (x[0] - _car2cyl.getCartesianOrigin()[0]);
1271 relPosition[1] = (x[1] - _car2cyl.getCartesianOrigin()[1]);
1272 relPosition[2] = (x[2] - _car2cyl.getCartesianOrigin()[2]);
1273
1274 T tmp[3] = {T(), T(), T()};
1275 _car2cyl(tmp, x);
1276 T rad = tmp[0];
1277 T phi = tmp[1];
1278
1279 T test = relPosition * normalAxis;
1280
1281 if ( (rad > _radWire && rad <= _cutoff * _radWire) &&
1282 (T(0) <= test && test <= _length) ) {
1283
1284 magneticFieldPolar[0] = _factor / util::pow(rad, 2)
1285 * (util::cos(phi));
1286 magneticFieldPolar[1] = _factor / util::pow(rad, 2) * util::sin(phi);
1287
1288 // changes radial and angular force components to cartesian components.
1289 outputTmp[0] = magneticFieldPolar[0] * util::cos(phi)
1290 - magneticFieldPolar[1] * util::sin(phi);
1291 outputTmp[1] = magneticFieldPolar[0] * util::sin(phi)
1292 + magneticFieldPolar[1] * util::cos(phi);
1293
1294 // if not in standard axis direction
1295 if ( !( util::nearZero(normalAxis[0]) && util::nearZero(normalAxis[1]) && util::nearZero(normalAxis[2] - 1) ) ) {
1296 Vector<T, 3> e3(T(), T(), T(1));
1297 Vector<T, 3> orientation =
1298 crossProduct3D(Vector<T, 3>(_car2cyl.getAxisDirection()), e3);
1299
1300 AngleBetweenVectors3D<T, S> angle(util::fromVector3(e3), util::fromVector3(orientation));
1301 T alpha[1] = {T()};
1302 T tmp2[3] = {T(), T(), T()};
1303 for (int i = 0; i < 3; ++i) {
1304 tmp2[i] = _car2cyl.getAxisDirection()[i];
1305 }
1306 angle(alpha, tmp2);
1307
1308 std::vector<T> origin(3, T());
1309 RotationRoundAxis3D<T, S> rotRAxis(origin, util::fromVector3(orientation), alpha[0]);
1310 rotRAxis(output, outputTmp);
1311 }
1312 else {
1313 output[0] = outputTmp[0];
1314 output[1] = outputTmp[1];
1315 output[2] = T();
1316 }
1317 }
1318 else {
1319 output[0] = outputTmp[0];
1320 output[1] = outputTmp[1];
1321 output[2] = outputTmp[1];
1322 }
1323
1324 return true;
1325}
1326
1327
1328} // end namespace olb
1329
1330#endif
M_PI
#define M_PI
Definition analyticalF.hh:43
olb::AnalyticalF
AnalyticalF are applications from DD to XD, where X is set by the constructor.
Definition analyticalBaseF.h:46
olb::AnalyticalPorousVelocity3D::AnalyticalPorousVelocity3D
AnalyticalPorousVelocity3D(SuperGeometry< T, 3 > &superGeometry, int material, T K_, T mu_, T gradP_, T radius_, T eps_=T(1))
Definition frameChangeF3D.hh:708
olb::AnalyticalPorousVelocity3D::operator()
bool operator()(T output[], const T input[]) override
Definition frameChangeF3D.hh:733
olb::AngleBetweenVectors3D
This class calculates the angle alpha between vector _referenceVector and any other vector x.
Definition frameChangeF3D.h:335
olb::AngleBetweenVectors3D::AngleBetweenVectors3D
AngleBetweenVectors3D(std::vector< T > referenceVector, std::vector< T > orientation)
constructor defines _referenceVector and _orientation
Definition frameChangeF3D.hh:760
olb::AngleBetweenVectors3D::operator()
bool operator()(T output[], const S x[]) override
operator writes angle between x and _referenceVector inro output field
Definition frameChangeF3D.hh:771
olb::CartesianToCylinder3D
This class converts Cartesian coordinates of point x to cylindrical coordinates wrote into output fie...
Definition frameChangeF3D.h:404
olb::CartesianToCylinder3D::getAxisDirection
olb::Vector< T, 3 > const & getAxisDirection() const
Read only access to _axisDirection.
Definition frameChangeF3D.hh:1020
olb::CartesianToCylinder3D::_cartesianOrigin
olb::Vector< T, 3 > _cartesianOrigin
origin of the Cartesian coordinate system
Definition frameChangeF3D.h:407
olb::CartesianToCylinder3D::getCartesianOrigin
olb::Vector< T, 3 > const & getCartesianOrigin() const
Read only access to _catesianOrigin.
Definition frameChangeF3D.hh:1034
olb::CartesianToCylinder3D::operator()
bool operator()(T output[], const S x[]) override
operator writes cylindrical coordinates of Cartesian point x into output field, output[0] = radius ( ...
Definition frameChangeF3D.hh:971
olb::CartesianToCylinder3D::CartesianToCylinder3D
CartesianToCylinder3D(olb::Vector< T, 3 > cartesianOrigin, T &angle, olb::Vector< T, 3 > orientation={T(1), T(), T()})
Definition frameChangeF3D.hh:916
olb::CartesianToCylinder3D::_axisDirection
olb::Vector< T, 3 > _axisDirection
direction of the axis along which the cylindrical coordinates are calculated
Definition frameChangeF3D.h:409
olb::CartesianToCylinder3D::setAngle
void setAngle(T angle)
Read and write access to _axisDirection using angle.
Definition frameChangeF3D.hh:1048
olb::CartesianToCylinder3D::_orientation
olb::Vector< T, 3 > _orientation
direction to know orientation for math positive to obtain angle phi of Cartesian point x
Definition frameChangeF3D.h:412
olb::CartesianToPolar2D
This class converts Cartesian coordinates of point x to polar coordinates wrote into output field (ou...
Definition frameChangeF2D.h:146
olb::CartesianToSpherical3D::operator()
bool operator()(T output[], const S x[]) override
operator writes shperical coordinates of Cartesian point x into output field, output[0] = radius ( >=...
Definition frameChangeF3D.hh:1101
olb::CartesianToSpherical3D::CartesianToSpherical3D
CartesianToSpherical3D(std::vector< T > cartesianOrigin, std::vector< T > axisDirection)
Definition frameChangeF3D.hh:1089
olb::CirclePoiseuille3D
Velocity profile for util::round pipes and a laminar flow of a Newtonian fluid: u(r)=u_max*(1-(r/R)^2...
Definition frameChangeF3D.h:215
olb::CirclePoiseuille3D::CirclePoiseuille3D
CirclePoiseuille3D(std::vector< T > axisPoint, std::vector< T > axisDirection, T maxVelocity, T radius, T scale=T(1))
Definition frameChangeF3D.hh:317
olb::CirclePoiseuilleStrainRate3D::operator()
bool operator()(T output[], const T input[]) override
Definition frameChangeF3D.hh:369
olb::CirclePoiseuilleStrainRate3D::CirclePoiseuilleStrainRate3D
CirclePoiseuilleStrainRate3D(UnitConverter< T, DESCRIPTOR > const &converter, T ly)
Definition frameChangeF3D.hh:359
olb::CirclePowerLaw3D
This functor returns a quadratic Poiseuille profile for use with a pipe with util::round cross-sectio...
Definition frameChangeF3D.h:154
olb::CirclePowerLaw3D::_normal
std::vector< T > _normal
Definition frameChangeF3D.h:157
olb::CirclePowerLaw3D::operator()
bool operator()(T output[], const T x[]) override
Definition frameChangeF3D.hh:233
olb::CirclePowerLaw3D::_center
olb::Vector< T, 3 > _center
Definition frameChangeF3D.h:156
olb::CirclePowerLaw3D::CirclePowerLaw3D
CirclePowerLaw3D(olb::Vector< T, 3 > axisPoint, std::vector< T > axisDirection, T maxVelocity, T radius, T n, T scale=T(1))
Definition frameChangeF3D.hh:166
olb::CirclePowerLawTurbulent3D
Velocity profile for util::round pipes and turbulent flows: u(r)=u_max*(1-r/R)^(1/n) The exponent n c...
Definition frameChangeF3D.h:194
olb::CirclePowerLawTurbulent3D::CirclePowerLawTurbulent3D
CirclePowerLawTurbulent3D(std::vector< T > axisPoint_, std::vector< T > axisDirection, T maxVelocity, T radius, T n=7, T turbulenceIntensity=0.05, T scale=T(1))
Definition frameChangeF3D.hh:242
olb::CirclePowerLawTurbulent3D::operator()
bool operator()(T output[], const T x[]) override
Definition frameChangeF3D.hh:283
olb::CylinderToCartesian3D::operator()
bool operator()(T output[], const S x[]) override
operator writes Cartesian coordinates of cylinder coordinates x[0] = radius >= 0, x[1] = phi in [0,...
Definition frameChangeF3D.hh:905
olb::CylinderToCartesian3D::CylinderToCartesian3D
CylinderToCartesian3D(std::vector< T > cylinderOrigin)
constructor defines _cylinderOrigin
Definition frameChangeF3D.hh:894
olb::EllipticPoiseuille3D::EllipticPoiseuille3D
EllipticPoiseuille3D(std::vector< T > center, T a, T b, T maxVel)
Definition frameChangeF3D.hh:683
olb::EllipticPoiseuille3D::operator()
bool operator()(T output[], const T x[]) override
Definition frameChangeF3D.hh:689
olb::GenericF::getName
std::string & getName()
read and write access to name
Definition genericF.hh:51
olb::MagneticFieldFromCylinder3D::_factor
T _factor
factor = mu0*4/3.*PI*radParticle^3*_Mp*radWire^2/r^3
Definition frameChangeF3D.h:538
olb::MagneticFieldFromCylinder3D::MagneticFieldFromCylinder3D
MagneticFieldFromCylinder3D(CartesianToCylinder3D< T, S > &car2cyl, T length, T radWire, T cutoff, T Mw)
Definition frameChangeF3D.hh:1246
olb::MagneticFieldFromCylinder3D::operator()
bool operator()(T output[], const S x[]) override
operator writes the magnetic force in a point x util::round a cylindrical wire into output field
Definition frameChangeF3D.hh:1261
olb::MagneticFieldFromCylinder3D::_Mw
T _Mw
saturation magnetization wire, linear scaling factor
Definition frameChangeF3D.h:536
olb::MagneticForceFromCylinder3D::_factor
T _factor
factor = mu0*4/3.*PI*radParticle^3*_Mp*radWire^2/r^3
Definition frameChangeF3D.h:516
olb::MagneticForceFromCylinder3D::_radParticle
T _radParticle
particle radius, cubic scaling factor
Definition frameChangeF3D.h:512
olb::MagneticForceFromCylinder3D::_Mw
T _Mw
saturation magnetization wire, linear scaling factor
Definition frameChangeF3D.h:508
olb::MagneticForceFromCylinder3D::MagneticForceFromCylinder3D
MagneticForceFromCylinder3D(CartesianToCylinder3D< T, S > &car2cyl, T length, T radWire, T cutoff, T Mw, T Mp=T(1), T radParticle=T(1), T mu0=4 *3.14159265e-7)
Magnetic field that creates magnetization in wire has to be orthogonal to the wire.
Definition frameChangeF3D.hh:1162
olb::MagneticForceFromCylinder3D::operator()
bool operator()(T output[], const S x[]) override
operator writes the magnetic force in a point x util::round a cylindrical wire into output field
Definition frameChangeF3D.hh:1180
olb::MagneticForceFromCylinder3D::_Mp
T _Mp
magnetization particle, linear scaling factor
Definition frameChangeF3D.h:510
olb::MagneticForceFromCylinder3D::_mu0
T _mu0
permeability constant, linear scaling factor
Definition frameChangeF3D.h:514
olb::OstreamManager
class for marking output with some text
Definition ostreamManager.h:81
olb::PolarToCartesian2D
This class converts polar coordinates of point x (x[0] = radius, x[1] = phi) to Cartesian coordinates...
Definition frameChangeF2D.h:124
olb::RectanglePoiseuille3D::RectanglePoiseuille3D
RectanglePoiseuille3D(olb::Vector< T, 3 > &x0_, olb::Vector< T, 3 > &x1_, olb::Vector< T, 3 > &x2_, std::vector< T > &maxVelocity_)
Definition frameChangeF3D.hh:401
olb::RectanglePoiseuille3D::operator()
bool operator()(T output[], const T x[]) override
Definition frameChangeF3D.hh:460
olb::RectangleTrigonometricPoiseuille3D::RectangleTrigonometricPoiseuille3D
RectangleTrigonometricPoiseuille3D(std::vector< T > &x0_, std::vector< T > &x1_, std::vector< T > &x2_, std::vector< T > &maxVelocity, int numOfPolynoms_=1)
Definition frameChangeF3D.hh:518
olb::RectangleTrigonometricPoiseuille3D::operator()
bool operator()(T output[], const T x[]) override
Definition frameChangeF3D.hh:614
olb::RotatingLinear3D::RotatingLinear3D
RotatingLinear3D(std::vector< T > axisPoint_, std::vector< T > axisDirection_, T w_, T scale_=1)
Definition frameChangeF3D.hh:46
olb::RotatingLinear3D::operator()
bool operator()(T output[], const T x[]) override
Definition frameChangeF3D.hh:53
olb::RotatingLinearAnnulus3D::RotatingLinearAnnulus3D
RotatingLinearAnnulus3D(std::vector< T > axisPoint_, std::vector< T > axisDirection_, T w_, T ri_, T ro_, T scale_=1)
Definition frameChangeF3D.hh:64
olb::RotatingLinearAnnulus3D::operator()
bool operator()(T output[], const T x[])
Definition frameChangeF3D.hh:71
olb::RotatingQuadratic1D::RotatingQuadratic1D
RotatingQuadratic1D(std::vector< T > axisPoint_, std::vector< T > axisDirection_, T w_, T scale_=1, T additive_=0)
Definition frameChangeF3D.hh:138
olb::RotatingQuadratic1D::axisPoint
std::vector< T > axisPoint
Definition frameChangeF3D.h:120
olb::RotatingQuadratic1D::axisDirection
std::vector< T > axisDirection
Definition frameChangeF3D.h:121
olb::RotatingQuadratic1D::operator()
bool operator()(T output[], const T x[]) override
Definition frameChangeF3D.hh:152
olb::RotationRoundAxis3D
This class saves coordinates of the resulting point after rotation in the output field.
Definition frameChangeF3D.h:356
olb::RotationRoundAxis3D::RotationRoundAxis3D
RotationRoundAxis3D(olb::Vector< T, 3 > origin, olb::Vector< T, 3 > rotAxisDirection, T alpha)
constructor defines _rotAxisDirection, _alpha and _origin
Definition frameChangeF3D.hh:836
olb::RotationRoundAxis3D::operator()
bool operator()(T output[], const S x[]) override
operator writes coordinates of the resulting point after rotation of x by angle _alpha in math positi...
Definition frameChangeF3D.hh:849
olb::SphericalToCartesian3D::operator()
bool operator()(T output[], const S x[]) override
operator writes spherical coordinates of Cartesian coordinates of x (x[0] = radius,...
Definition frameChangeF3D.hh:1079
olb::SuperGeometry
Representation of a statistic for a parallel 2D geometry.
Definition superGeometry.h:70
olb::SuperGeometry::getStatistics
SuperGeometryStatistics< T, D > & getStatistics()
Returns the statistics object.
Definition superGeometry.hh:190
olb::UnitConverter
Conversion between physical and lattice units, as well as discretization.
Definition unitConverter.h:85
olb::UnitConverter::getCharPhysVelocity
constexpr T getCharPhysVelocity() const
return characteristic velocity in physical units
Definition unitConverter.h:148
olb::Vector
Plain old scalar vector.
Definition vector.h:47
frameChangeF2D.h
This file contains two different classes of functors, in the FIRST part.
frameChangeF3D.h
This file contains two different classes of functors, in the FIRST part.
olb::util::sqrt
cpu::simd::Pack< T > sqrt(cpu::simd::Pack< T > value)
Definition pack.h:100
olb::util::fromVector3
std::vector< T > fromVector3(const Vector< T, 3 > &vec)
Definition vectorHelpers.h:126
olb::util::atan2
ADf< T, DIM > atan2(const T &y, const ADf< T, DIM > &x)
Definition aDiff.h:623
olb::util::sin
ADf< T, DIM > sin(const ADf< T, DIM > &a)
Definition aDiff.h:569
olb::util::pow
cpu::simd::Pack< T > pow(cpu::simd::Pack< T > base, cpu::simd::Pack< T > exp)
Definition pack.h:112
olb::util::cosh
ADf< T, DIM > cosh(const ADf< T, DIM > &a)
Definition aDiff.h:657
olb::util::fabs
cpu::simd::Pack< T > fabs(cpu::simd::Pack< T > value)
Definition pack.h:106
olb::util::normalize
Vector< T, D > normalize(const Vector< T, D > &a)
Definition vectorHelpers.h:193
olb::util::acos
ADf< T, DIM > acos(const ADf< T, DIM > &a)
Definition aDiff.h:605
olb::util::nearZero
bool nearZero(const ADf< T, DIM > &a)
Definition aDiff.h:1087
olb::util::tanh
ADf< T, DIM > tanh(const ADf< T, DIM > &a)
Definition aDiff.h:663
olb::util::atan
ADf< T, DIM > atan(const ADf< T, DIM > &a)
Definition aDiff.h:614
olb::util::cos
ADf< T, DIM > cos(const ADf< T, DIM > &a)
Definition aDiff.h:578
olb
Top level namespace for all of OpenLB.
Definition boundaryPostProcessors2D.h:34
olb::norm
constexpr T norm(const ScalarVector< T, D, IMPL > &a)
Euclidean vector norm.
Definition scalarVector.h:66
olb::crossProduct3D
constexpr Vector< T, 3 > crossProduct3D(const ScalarVector< T, 3, IMPL > &a, const ScalarVector< T, 3, IMPL_ > &b) any_platform
Definition vector.h:224
olb::normalize
constexpr Vector< T, D > normalize(const ScalarVector< T, D, IMPL > &a, T scale=T{1})
Definition vector.h:245
superGeometry.h
Representation of a parallel 2D geometry – header file.
Generated on Tue Mar 5 2024 09:57:48 for OpenLB by doxygen 1.10.0