LifeV
structure/examples/example_evaluatingScalarVectorialTensorialQuantitiesUsingETA/ud_functions.cpp
Go to the documentation of this file.
1 //@HEADER
2 /*
3 *******************************************************************************
4 
5  Copyright (C) 2004, 2005, 2007 EPFL, Politecnico di Milano, INRIA
6  Copyright (C) 2010 EPFL, Politecnico di Milano, Emory University
7 
8  This file is part of LifeV.
9 
10  LifeV is free software; you can redistribute it and/or modify
11  it under the terms of the GNU Lesser General Public License as published by
12  the Free Software Foundation, either version 3 of the License, or
13  (at your option) any later version.
14 
15  LifeV is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  Lesser General Public License for more details.
19 
20  You should have received a copy of the GNU Lesser General Public License
21  along with LifeV. If not, see <http://www.gnu.org/licenses/>.
22 
23 *******************************************************************************
24 */
25 //@HEADER
26 
27 /*!
28  * @file
29  * @brief File containing the boundary conditions for the Monolithic Test
30  *
31  * @date 2009-04-09
32  * @author Paolo Tricerri <paolo.tricerri@epfl.ch>
33  *
34  * @maintainer Paolo Tricerri <paolo.tricerri@epfl.ch>
35  *
36  * Contains the functions to be assigned as boundary conditions, in the file boundaryConditions.hpp . The functions
37  * can depend on time and space, while they can take in input an ID specifying one of the three principal axis
38  * if the functions to assign is vectorial and the boundary condition is of type \c Full \c.
39  */
40 
41 #include "ud_functions.hpp"
42 #include<lifev/core/array/VectorSmall.hpp>
43 #include<lifev/core/array/MatrixSmall.hpp>
44 
45 #define PI 3.14159265359
46 
47 namespace LifeV
48 {
49 
50 Real f (const Real& /*t*/, const Real& /*x*/, const Real& /*y*/, const Real& /*z*/, const ID& i)
51 {
52  switch (i)
53  {
54  case 0:
55  return 0.0;
56  break;
57  case 1:
58  return 0.0;
59  break;
60  case 2:
61  return 0.0;
62  break;
63  default:
64  ERROR_MSG ("This entrie is not allowed: ud_functions.hpp");
65  return 0.;
66  break;
67  }
68 }
69 
70 
71 Real fzero_scalar (const Real& /*t*/, const Real& /*x*/, const Real& /*y*/, const Real& /*z*/, const ID& /*i*/)
72 {
73  return 0.0;
74 }
75 
76 Real InternalPressure (const Real& /*t*/, const Real& /*x*/, const Real& /*y*/, const Real& /*z*/, const ID& /*i*/)
77 {
78  return -1e+5;
79  //return -260000*sin(80*3.141592*t);
80 }
81 
82 // Initial displacement and velocity
83 Real d0 (const Real& /*t*/, const Real& x, const Real& y, const Real& z, const ID& i)
84 {
85 
86  switch (i)
87  {
88  case 0:
89  return 0.088002 * ( x + 0.5 );
90  break;
91  case 1:
92  return - ( 0.02068 * 2.0 ) * ( y );
93  break;
94  case 2:
95  return - ( 0.02068 * 2.0 ) * ( z );
96  break;
97  default:
98  ERROR_MSG ("This entry is not allowed: ud_functions.hpp");
99  return 0.;
100  break;
101  }
102 
103 }
104 
105 Real w0 (const Real& /*t*/, const Real& /*x*/, const Real& /*y*/, const Real& /*z*/, const ID& i)
106 {
107 
108  switch (i)
109  {
110  case 0:
111  return 0.0;
112  break;
113  case 1:
114  return 0.0;
115  break;
116  case 2:
117  return 0.0;
118  break;
119  default:
120  ERROR_MSG ("This entrie is not allowed: ud_functions.hpp");
121  return 0.;
122  break;
123  }
124 }
125 
126 Real a0 (const Real& /*t*/, const Real& /*x*/, const Real& /*y*/, const Real& /*z*/, const ID& i)
127 {
128 
129  switch (i)
130  {
131  case 0:
132  return 0.0;
133  break;
134  case 1:
135  return 0.0;
136  break;
137  case 2:
138  return 0.0;
139  break;
140  default:
141  ERROR_MSG ("This entrie is not allowed: ud_functions.hpp");
142  return 0.;
143  break;
144  }
145 }
146 
147 
148 //----------------------------------------------Boundary Conditions--------------
149 Real bcZero (const Real& /*t*/, const Real& /*X*/, const Real& /*Y*/, const Real& /*Z*/, const ID& /*i*/)
150 {
151  return 0.;
152 }
153 
154 Real bcNonZero (const Real& t, const Real& X, const Real& Y, const Real& Z, const ID& i)
155 {
156  //Real pressure(200000);
157 
158  // Real top = 1000000;
159 
160  // return top * Y;
161 
162  return 10000;
163 
164  // Real highestPressure(6.666e+6);
165  // Real totalTime = 20.0;
166  // Real halfTime = totalTime / 2.0;
167 
168  // Real a = ( highestPressure / 2 ) * ( 1/ ((totalTime/2)*(totalTime/2)) );
169 
170  // if ( t <= halfTime )
171  // pressure = a * t*t;
172 
173  // if ( t > halfTime )
174  // pressure = - a * (t - totalTime)*(t - totalTime) + highestPressure;
175 
176  // switch (i)
177  // {
178  // case 0:
179  // return 0.0;
180  // break;
181  // case 1:
182  // return pressure;
183  // break;
184  // case 2:
185  // return 0.0;
186  // break;
187  // default:
188  // ERROR_MSG ("This entrie is not allowed: ud_functions.hpp");
189  // return 0.;
190  // break;
191  // }
192 }
193  Real smoothPressure(const Real& t, const Real& x, const Real& y, const Real& /*Z*/, const ID& i)
194  {
195  Real radius = std::sqrt( x*x + y*y);
196  Real pressure(0);
197 
198  Real highestPressure(200000);
199  Real totalTime = 4.5;
200  Real halfTime = totalTime / 2.0;
201 
202  Real a = ( highestPressure / 2 ) * ( 1/ ( halfTime*halfTime ) );
203 
204  if ( t <= halfTime )
205  pressure = a * t*t;
206 
207  if ( t > halfTime )
208  pressure = - a * (t - totalTime)*(t - totalTime) + highestPressure;
209 
210  switch (i)
211  {
212  case 0:
213  return pressure * ( x / radius ) ;
214  break;
215  case 1:
216  return pressure * ( y / radius ) ;
217  break;
218  case 2:
219  return 0.0;
220  break;
221 
222  }
223  return 0;
224 
225  }
226 
227 Real traction (const Real& /*t*/, const Real& /*X*/, const Real& /*Y*/, const Real& /*Z*/, const ID& /*i*/)
228 {
229  return 10000.;
230 }
231 
232 Real referenceDirection ( const Real& /*t*/, const Real& x, const Real& y, const Real& z, const ID& i)
233 {
234  Real theta = 0.0; // value for anisotropic characterization taken from Robertson // ( PI / 6.0 );
235  Real thetaChangeOfVariable = std::atan( y / x );
236 
237  if( x < 0 )
238  {
239  // This is due to the periodicity of std::atan ( ref. official documentation )
240  thetaChangeOfVariable += PI;
241  }
242 
243  // Half thorus
244  // Real thetaFiber ( theta );
245  // Real xT (0.138);
246  // Real rT ( xT );
247  // Real thetaFree(0);
248  // Real thetaRotation(0);
249  // Real thetaPosition(0);
250  // thetaFree = std::atan( z / ( x - xT ) );
251  // thetaPosition = thetaFree;
252  // if ( x < xT )
253  // {
254  // thetaPosition = thetaFree + PI;
255  // thetaRotation = PI/2.0 - std::fabs( thetaFree );
256  // }
257  // else
258  // {
259  // thetaRotation = thetaFree - PI/2.0;
260  // }
261 
262  // Real xCenter;
263  // Real zCenter;
264  // xCenter = rT * std::cos( thetaPosition ) + xT;
265  // zCenter = rT * std::sin( thetaPosition );
266 
267  // VectorSmall<3> positionCenter;
268  // positionCenter[0] = xCenter;
269  // positionCenter[1] = 0;
270  // positionCenter[2] = zCenter;
271 
272  // VectorSmall<3> position;
273  // position[0] = x;
274  // position[1] = y;
275  // position[2] = z;
276 
277  // MatrixSmall<3,3> changeOfVariable;
278  // changeOfVariable(0,0) = std::cos( thetaRotation ); changeOfVariable(0,1) = 0.0; changeOfVariable(0,2) = std::sin( thetaRotation );
279  // changeOfVariable(1,0) = 0.0; changeOfVariable(1,1) = 1.0; changeOfVariable(1,2) = 0.0;
280  // changeOfVariable(2,0) = -std::sin( thetaRotation ); changeOfVariable(2,1) = 0.0; changeOfVariable(2,2) = std::cos( thetaRotation );
281 
282  // VectorSmall<3> localPosition;
283  // localPosition = changeOfVariable * ( position - positionCenter );
284 
285  // Real thetaPositionOnSection(0);
286  // thetaPositionOnSection = std::atan( localPosition[2] / localPosition[1] );
287 
288  // if ( localPosition[1] < 0 )
289  // thetaPositionOnSection += PI;
290 
291  // VectorSmall<3> localFibers;
292  // localFibers[0] = std::sin( thetaFiber );
293  // localFibers[1] = std::cos( thetaFiber ) * std::sin( thetaPositionOnSection );
294  // localFibers[2] = - std::cos( thetaFiber ) * std::cos( thetaPositionOnSection );
295 
296  // VectorSmall<3> originFibers;
297  // originFibers = changeOfVariable.transpose() * localFibers + positionCenter;
298 
299  // originFibers.normalize();
300 
301  switch (i)
302  {
303  case 0:
304  // Tube
305  return - std::sin( thetaChangeOfVariable ) * std::cos( theta );
306  // half thorus
307  //return originFibers[0];
308  // Cube
309  // return std::sin( theta );
310  break;
311  case 1:
312  // Tube
313  return std::cos( thetaChangeOfVariable ) * std::cos( theta );
314  // half thorus
315  //return originFibers[1];
316  // Cube
317  //return std::cos( theta );
318  break;
319  case 2:
320  // Tube
321  return std::sin( theta );
322  // half thorus
323  //return originFibers[2];
324  // Cube
325  // return 0.0;
326  break;
327  default:
328  ERROR_MSG ("This entrie is not allowed: ud_functions.hpp");
329  return 0.;
330  break;
331  }
332 }
333 
334 
335 
336 //----------------------------------------------Fibers Directions--------------
337 Real Family1 ( const Real& /*t*/, const Real& x, const Real& y, const Real& z, const ID& i)
338 {
339  Real theta = 0.74725; // value for anisotropic characterization taken from Robertson // ( PI / 6.0 );
340  Real thetaChangeOfVariable = std::atan( y / x );
341 
342  if( x < 0 )
343  {
344  // This is due to the periodicity of std::atan ( ref. official documentation )
345  thetaChangeOfVariable += PI;
346  }
347 
348  // Half thorus
349  // Real thetaFiber ( theta );
350  // Real xT (0.138);
351  // Real rT ( xT );
352  // Real thetaFree(0);
353  // Real thetaRotation(0);
354  // Real thetaPosition(0);
355  // thetaFree = std::atan( z / ( x - xT ) );
356  // thetaPosition = thetaFree;
357  // if ( x < xT )
358  // {
359  // thetaPosition = thetaFree + PI;
360  // thetaRotation = PI/2.0 - std::fabs( thetaFree );
361  // }
362  // else
363  // {
364  // thetaRotation = thetaFree - PI/2.0;
365  // }
366 
367  // Real xCenter;
368  // Real zCenter;
369  // xCenter = rT * std::cos( thetaPosition ) + xT;
370  // zCenter = rT * std::sin( thetaPosition );
371 
372  // VectorSmall<3> positionCenter;
373  // positionCenter[0] = xCenter;
374  // positionCenter[1] = 0;
375  // positionCenter[2] = zCenter;
376 
377  // VectorSmall<3> position;
378  // position[0] = x;
379  // position[1] = y;
380  // position[2] = z;
381 
382  // MatrixSmall<3,3> changeOfVariable;
383  // changeOfVariable(0,0) = std::cos( thetaRotation ); changeOfVariable(0,1) = 0.0; changeOfVariable(0,2) = std::sin( thetaRotation );
384  // changeOfVariable(1,0) = 0.0; changeOfVariable(1,1) = 1.0; changeOfVariable(1,2) = 0.0;
385  // changeOfVariable(2,0) = -std::sin( thetaRotation ); changeOfVariable(2,1) = 0.0; changeOfVariable(2,2) = std::cos( thetaRotation );
386 
387  // VectorSmall<3> localPosition;
388  // localPosition = changeOfVariable * ( position - positionCenter );
389 
390  // Real thetaPositionOnSection(0);
391  // thetaPositionOnSection = std::atan( localPosition[2] / localPosition[1] );
392 
393  // if ( localPosition[1] < 0 )
394  // thetaPositionOnSection += PI;
395 
396  // VectorSmall<3> localFibers;
397  // localFibers[0] = std::sin( thetaFiber );
398  // localFibers[1] = std::cos( thetaFiber ) * std::sin( thetaPositionOnSection );
399  // localFibers[2] = - std::cos( thetaFiber ) * std::cos( thetaPositionOnSection );
400 
401  // VectorSmall<3> originFibers;
402  // originFibers = changeOfVariable.transpose() * localFibers + positionCenter;
403 
404  // originFibers.normalize();
405 
406  switch (i)
407  {
408  case 0:
409  // Tube
410  return - std::sin( thetaChangeOfVariable ) * std::cos( theta );
411  // half thorus
412  //return originFibers[0];
413  // Cube
414  // return std::sin( theta );
415  break;
416  case 1:
417  // Tube
418  return std::cos( thetaChangeOfVariable ) * std::cos( theta );
419  // half thorus
420  //return originFibers[1];
421  // Cube
422  //return std::cos( theta );
423  break;
424  case 2:
425  // Tube
426  return std::sin( theta );
427  // half thorus
428  //return originFibers[2];
429  // Cube
430  // return 0.0;
431  break;
432  default:
433  ERROR_MSG ("This entrie is not allowed: ud_functions.hpp");
434  return 0.;
435  break;
436  }
437 }
438 
439 Real Family2 ( const Real& /*t*/, const Real& x, const Real& y, const Real& z, const ID& i)
440 {
441  Real theta = -0.74725;
442  Real thetaChangeOfVariable = std::atan( y / x );
443 
444  if( x < 0 )
445  {
446  // This is due to the periodicity of std::atan ( ref. official documentation )
447  thetaChangeOfVariable += PI;
448  }
449 
450  // Half thorus
451  // Real thetaFiber ( theta );
452  // Real xT (0.138);
453  // Real rT ( xT );
454  // Real thetaFree(0);
455  // Real thetaRotation(0);
456  // Real thetaPosition(0);
457  // thetaFree = std::atan( z / ( x - xT ) );
458  // thetaPosition = thetaFree;
459  // if ( x < xT )
460  // {
461  // thetaPosition = thetaFree + PI;
462  // thetaRotation = PI/2.0 - std::fabs( thetaFree );
463  // }
464  // else
465  // {
466  // thetaRotation = thetaFree - PI/2.0;
467  // }
468 
469  // Real xCenter;
470  // Real zCenter;
471  // xCenter = rT * std::cos( thetaPosition ) + xT;
472  // zCenter = rT * std::sin( thetaPosition );
473 
474  // VectorSmall<3> positionCenter;
475  // positionCenter[0] = xCenter;
476  // positionCenter[1] = 0;
477  // positionCenter[2] = zCenter;
478 
479  // VectorSmall<3> position;
480  // position[0] = x;
481  // position[1] = y;
482  // position[2] = z;
483 
484  // MatrixSmall<3,3> changeOfVariable;
485  // changeOfVariable(0,0) = std::cos( thetaRotation ); changeOfVariable(0,1) = 0.0; changeOfVariable(0,2) = std::sin( thetaRotation );
486  // changeOfVariable(1,0) = 0.0; changeOfVariable(1,1) = 1.0; changeOfVariable(1,2) = 0.0;
487  // changeOfVariable(2,0) = -std::sin( thetaRotation ); changeOfVariable(2,1) = 0.0; changeOfVariable(2,2) = std::cos( thetaRotation );
488 
489  // VectorSmall<3> localPosition;
490  // localPosition = changeOfVariable * ( position - positionCenter );
491 
492  // Real thetaPositionOnSection(0);
493  // thetaPositionOnSection = std::atan( localPosition[2] / localPosition[1] );
494 
495  // if ( localPosition[1] < 0 )
496  // thetaPositionOnSection += PI;
497 
498  // VectorSmall<3> localFibers;
499  // localFibers[0] = std::sin( thetaFiber );
500  // localFibers[1] = std::cos( thetaFiber ) * std::sin( thetaPositionOnSection );
501  // localFibers[2] = - std::cos( thetaFiber ) * std::cos( thetaPositionOnSection );
502 
503  // VectorSmall<3> originFibers;
504  // originFibers = changeOfVariable.transpose() * localFibers + positionCenter;
505 
506  // originFibers.normalize();
507 
508  switch (i)
509  {
510  case 0:
511  // Tube
512  return - std::sin( thetaChangeOfVariable ) * std::cos( theta );
513  // half thorus
514  //return originFibers[0];
515  // Cube
516  //return std::sin( theta );
517  break;
518  case 1:
519  // Tube
520  return std::cos( thetaChangeOfVariable ) * std::cos( theta );
521  // half thorus
522  //return originFibers[1];
523  // Cube
524  //return std::cos( theta );
525  break;
526  case 2:
527  // Tube
528  return std::sin( theta );
529  // half thorus
530  //return originFibers[2];
531  // Cube
532  //return 0.0;
533  break;
534  default:
535  ERROR_MSG ("This entrie is not allowed: ud_functions.hpp");
536  return 0.;
537  break;
538  }
539 }
540 
541 Real Family3 ( const Real& /*t*/, const Real& x, const Real& y, const Real& z, const ID& i)
542 {
543 
544  switch (i)
545  {
546  case 0:
547  return 0.0;
548  break;
549  case 1:
550  return 0.0;
551  break;
552  case 2:
553  return -1.0;
554  break;
555  default:
556  ERROR_MSG ("This entrie is not allowed: ud_functions.hpp");
557  return 0.;
558  break;
559  }
560 }
561 
562 
563 Real Family4 ( const Real& /*t*/, const Real& x, const Real& y, const Real& z, const ID& i)
564 {
565 
566  switch (i)
567  {
568  case 0:
569  return 0.0;
570  break;
571  case 1:
572  return 0.0;
573  break;
574  case 2:
575  return 0.0;
576  break;
577  default:
578  ERROR_MSG ("This entrie is not allowed: ud_functions.hpp");
579  return 0.;
580  break;
581  }
582 }
583 
584 
585 Real Family5 ( const Real& /*t*/, const Real& x, const Real& y, const Real& z, const ID& i)
586 {
587 
588  switch (i)
589  {
590  case 0:
591  return 0.0;
592  break;
593  case 1:
594  return 0.0;
595  break;
596  case 2:
597  return 0.0;
598  break;
599  default:
600  ERROR_MSG ("This entrie is not allowed: ud_functions.hpp");
601  return 0.;
602  break;
603  }
604 }
605 
606 Real Family6 ( const Real& /*t*/, const Real& x, const Real& y, const Real& z, const ID& i)
607 {
608 
609  switch (i)
610  {
611  case 0:
612  return 0.0;
613  break;
614  case 1:
615  return 0.0;
616  break;
617  case 2:
618  return 0.0;
619  break;
620  default:
621  ERROR_MSG ("This entrie is not allowed: ud_functions.hpp");
622  return 0.;
623  break;
624  }
625 }
626 
627 // Method for the definition of the fibers
628 fibersDirectionList::fibersDirectionList() :
629  M_mapNameDefinition( )
630 {}
631 
632 fibersDirectionList::~fibersDirectionList()
633 {}
634 
635 void fibersDirectionList::setupFiberDefinitions( const UInt nbFamilies )
636 {
637  // At the moment the creation of the table of fiber functions is done
638  // manually. There should be a way to make it automatically. Btw, only
639  // the first nbFamilies that are set in the data file are taken into account
640 
641  ASSERT( nbFamilies < 6, "At the moment, a maximum number = 6 of families can be used! If you want more \n modifiy the file ud_functions.hpp in the application folder." );
642 
643  // Creation of the database of functions
644  fiberFunctionPtr_Type pointerToFunction( new fiberFunction_Type( Family1 ) );
645  M_mapNameDefinition.insert( std::pair<std::string, fiberFunctionPtr_Type>
646  ( "Family1", pointerToFunction ) );
647 
648  pointerToFunction.reset( new fiberFunction_Type( Family2 ) );
649  M_mapNameDefinition.insert( std::pair<std::string, fiberFunctionPtr_Type>
650  ( "Family2", pointerToFunction ) );
651 
652  pointerToFunction.reset( new fiberFunction_Type( Family3 ) );
653  M_mapNameDefinition.insert( std::pair<std::string, fiberFunctionPtr_Type>
654  ( "Family3", pointerToFunction ) );
655 
656  pointerToFunction.reset( new fiberFunction_Type( Family4 ) );
657  M_mapNameDefinition.insert( std::pair<std::string, fiberFunctionPtr_Type>
658  ( "Family4", pointerToFunction ) );
659 
660  pointerToFunction.reset( new fiberFunction_Type( Family5 ) );
661  M_mapNameDefinition.insert( std::pair<std::string, fiberFunctionPtr_Type>
662  ( "Family5", pointerToFunction ) );
663 
664  pointerToFunction.reset( new fiberFunction_Type( Family6 ) );
665  M_mapNameDefinition.insert( std::pair<std::string, fiberFunctionPtr_Type>
666  ( "Family6", pointerToFunction ) );
667 
668 
669 }
670 
671 fibersDirectionList::fiberFunctionPtr_Type fibersDirectionList::fiberDefinition( const std::string nameFamily )
672 {
673 
674  mapNameDefinitionFiberFunction_Type::const_iterator IT;
675 
676  IT = M_mapNameDefinition.find ( nameFamily );
677 
678  if ( IT != M_mapNameDefinition.end() )
679  {
680  return IT->second;
681  }
682  else
683  {
684  std::cout << " Wrong identification of the fiber function! " << std::endl;
685  fiberFunctionPtr_Type pointerToFunction( new fiberFunction_Type() );
686 
687  return pointerToFunction;
688  }
689 }
690 
691 }
Real d0(const Real &t, const Real &x, const Real &y, const Real &z, const ID &i)
#define PI
void updateInverseJacobian(const UInt &iQuadPt)
Real w0(const Real &t, const Real &x, const Real &y, const Real &z, const ID &i)
#define ERROR_MSG(A)
Definition: LifeAssert.hpp:69
#define ASSERT(X, A)
Definition: LifeAssert.hpp:90
uint32_type ID
IDs.
Definition: LifeV.hpp:194
Real f(const Real &t, const Real &x, const Real &y, const Real &z, const ID &i)
double Real
Generic real data.
Definition: LifeV.hpp:175
uint32_type UInt
generic unsigned integer (used mainly for addressing)
Definition: LifeV.hpp:191