LifeV
OneDFSIData.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 a class for 1D model data handling.
30  *
31  * @version 1.0
32  * @date 01-07-2004
33  * @author Vincent Martin
34  *
35  * @version 2.0
36  * @date 12-04-2010
37  * @author Cristiano Malossi <cristiano.malossi@epfl.ch>
38  *
39  * @contributors Ricardo Ruiz-Baier <ricardo.ruiz@epfl.ch>
40  * @maintainer Cristiano Malossi <cristiano.malossi@epfl.ch>
41  */
42 
43 #include <lifev/one_d_fsi/solver/OneDFSIData.hpp>
44 
45 namespace LifeV
46 {
47 
48 // ===================================================
49 // Constructors
50 // ===================================================
51 OneDFSIData::OneDFSIData() :
52  M_physicsType (),
53  M_fluxType (),
54  M_sourceType (),
55  M_timeDataPtr (),
56  M_meshPtr ( new mesh_Type ),
57  M_viscoelasticWall (),
58  M_viscoelasticAngle (),
59  M_viscoelasticPeriod (),
60  M_viscoelasticCoefficient (),
61  M_inertialWall (),
62  M_densityWall (),
63  M_inertialModulus (),
64  M_longitudinalWall (),
65  M_postprocessingDirectory (),
66  M_postprocessingFile (),
67  M_CFLmax (),
68  M_jacobianPerturbationArea (),
69  M_jacobianPerturbationFlowRate(),
70  M_jacobianPerturbationStress(),
71  M_computeCoefficients (),
72  M_powerLawCoefficient (),
73  M_density (),
74  M_viscosity (),
75  M_thickVessel (),
76  M_young (),
77  M_poisson (),
78  M_externalPressure (),
79  M_venousPressure (),
80  M_robertsonCorrection (),
81  M_thickness (),
82  M_friction (),
83  M_area0 (),
84  M_alpha (),
85  M_beta0 (),
86  M_beta1 (),
87  M_dArea0dz (),
88  M_dAlphadz (),
89  M_dBeta0dz (),
90  M_dBeta1dz (),
91  M_flux11 (),
92  M_flux12 (),
93  M_flux21 (),
94  M_flux22 (),
95  M_celerity1 (),
96  M_celerity2 (),
97  M_celerity1LeftEigenvector1 (),
98  M_celerity1LeftEigenvector2 (),
99  M_celerity2LeftEigenvector1 (),
100  M_celerity2LeftEigenvector2 (),
101  M_source10 (),
102  M_source20 (),
103  M_source11 (),
104  M_source12 (),
105  M_source21 (),
106  M_source22 ()
107 {
108 }
109 
110 // ===================================================
111 // Methods
112 // ===================================================
113 void
114 OneDFSIData::setup ( const GetPot& dataFile, const std::string& section )
115 {
116  // Model Type
117  M_physicsType = OneDFSI::physicsMap[ dataFile ( ( section + "/Model/PhysicsType" ).data(), "OneD_1DLinearPhysics" ) ];
118  M_fluxType = OneDFSI::fluxMap[ dataFile ( ( section + "/Model/FluxType" ).data(), "OneD_1DLinearFlux" ) ];
119  M_sourceType = OneDFSI::sourceMap[ dataFile ( ( section + "/Model/SourceType" ).data(), "OneD_1DLinearSource" ) ];
120 
121  // If data time has not been set
122  if ( !M_timeDataPtr.get() )
123  {
124  M_timeDataPtr.reset ( new time_Type ( dataFile, (section + "/time_discretization" ).data() ) );
125  }
126 
127  // Mesh setup - Space Discretization
128  Real length = dataFile ( ( section + "/space_discretization/Length" ).data(), 1. );
129  Real numberOfElements = dataFile ( ( section + "/space_discretization/NumberOfElements" ).data(), 10 );
130 
131  regularMesh1D ( *M_meshPtr, 1, numberOfElements, false, length, 0 );
132 
133  //std::cout << " 1D- Mesh nodes: " << M_meshPtr->numPoints() << std::endl;
134  //std::cout << " 1D- Mesh elements: " << M_meshPtr->numElements() << std::endl;
135 
136  // Physical Wall
137  M_viscoelasticWall = dataFile ( ( section + "/PhysicalWall/ViscoelasticWall" ).data(), false );
138  M_viscoelasticAngle = dataFile ( ( section + "/PhysicalWall/ViscoelasticAngle" ).data(), 5. ) * M_PI / 180.;
139  M_viscoelasticPeriod = dataFile ( ( section + "/PhysicalWall/ViscoelasticPeriod" ).data(), 0.3 * 0.8 );
140 
141  M_inertialWall = dataFile ( ( section + "/PhysicalWall/InertialWall" ).data(), false );
142  M_densityWall = dataFile ( ( section + "/PhysicalWall/DensityWall" ).data(), 1. );
143  M_inertialModulus = dataFile ( ( section + "/PhysicalWall/coeffA" ).data(), 0. );
144 
145  M_longitudinalWall = dataFile ( ( section + "/PhysicalWall/LongitudinalWall" ).data(), false );
146 
147  // Miscellaneous
148  M_postprocessingDirectory = dataFile ( ( section + "/miscellaneous/post_dir" ).data(), "./" );
149  M_postprocessingFile = dataFile ( ( section + "/miscellaneous/post_file" ).data(), "sol" );
150  M_CFLmax = dataFile ( ( section + "/miscellaneous/CFLmax" ).data(), std::sqrt (3.) / 3. );
151 
152  if ( M_CFLmax > std::sqrt (3.) / 3. )
153  {
154  std::cout << "!!! WARNING: CFLmax greater than the theoretical value (see MOX21, eq. 1.47) - CONVERGENCE NOT GUARANTEED !!!" << std::endl;
155  }
156 
157  // Jacobian perturbation
158  M_jacobianPerturbationArea = dataFile ( ( section + "/JacobianPerturbation/deltaArea" ).data(), 0.001 );
159  M_jacobianPerturbationFlowRate = dataFile ( ( section + "/JacobianPerturbation/deltaFlowRate" ).data(), 0.001 );
160  M_jacobianPerturbationStress = dataFile ( ( section + "/JacobianPerturbation/deltaStress" ).data(), 1.0 );
161 
162  // Physical Parameters
163  M_computeCoefficients = dataFile ( ( section + "/PhysicalParameters/ComputeCoefficients" ).data(), false );
164  M_powerLawCoefficient = dataFile ( ( section + "/PhysicalParameters/PowerlawCoefficient" ).data(), 2 );
165 
166  M_density = dataFile ( ( section + "/PhysicalParameters/density" ).data(), 1. );
167  M_viscosity = dataFile ( ( section + "/PhysicalParameters/viscosity" ).data(), 0.035 );
168 
169  M_thickVessel = dataFile ( ( section + "/PhysicalParameters/thickVessel" ).data(), false );
170  M_young = dataFile ( ( section + "/PhysicalParameters/young" ).data(), 4.0E6 );
171  M_poisson = dataFile ( ( section + "/PhysicalParameters/poisson" ).data(), 0.5 );
172 
173  M_externalPressure = dataFile ( ( section + "/PhysicalParameters/externalPressure" ).data(), 0. );
174  M_venousPressure = dataFile ( ( section + "/PhysicalParameters/venousPressure" ).data(), 0. );
175  M_friction = dataFile ( ( section + "/PhysicalParameters/Kr" ).data(), 1. );
176 
177  M_robertsonCorrection = dataFile ( ( section + "/PhysicalParameters/RobertsonCorrection" ).data(), 1. );
178 
179  M_computeCoefficients = dataFile ( ( section + "/PhysicalParameters/ComputeCoefficients" ).data(), false );
180 
181  // Reset all the containers
182  resetContainers();
183 
184  // Select a law for the coefficients
185  std::map< std::string, OneD_distributionLaw > distributionLawMap;
186  distributionLawMap["uniform"] = uniform;
187  distributionLawMap["linear"] = linear;
188  distributionLawMap["pointwise"] = pointwise;
189 
190  OneD_distributionLaw distributionLaw = distributionLawMap[ dataFile ( ( section + "/PhysicalParameters/DistributionLaw" ).data(), "uniform" ) ];
191  switch ( distributionLaw )
192  {
193  case uniform:
194 
195  for ( UInt i = 0; i < M_meshPtr->numPoints() ; ++i )
196  {
197  // Physical Parameters
198  M_thickness[i] = dataFile ( ( section + "/PhysicalParameters/thickness" ).data(), 0. );
199  M_viscoelasticCoefficient[i] = dataFile ( ( section + "/PhysicalWall/ViscoelasticCoefficient" ).data(), 0. );
200 
201  M_area0[i] = dataFile ( ( section + "/PhysicalParameters/Area0" ).data(), M_PI );
202  M_alpha[i] = dataFile ( ( section + "/PhysicalParameters/AlphaCoriolis" ).data(), 1. / M_robertsonCorrection );
203  M_beta0[i] = dataFile ( ( section + "/PhysicalParameters/Beta0" ).data(), 1.e6 );
204  M_beta1[i] = dataFile ( ( section + "/PhysicalParameters/Beta1" ).data(), 0.5 );
205 
206  // Linear Parameters
207  M_flux11[i] = dataFile ( ( section + "/LinearParameters/Flux11" ).data(), 1. );
208  M_flux12[i] = dataFile ( ( section + "/LinearParameters/Flux12" ).data(), 0. );
209  M_flux21[i] = dataFile ( ( section + "/LinearParameters/Flux21" ).data(), 0. );
210  M_flux22[i] = dataFile ( ( section + "/LinearParameters/Flux22" ).data(), 1. );
211  M_celerity1[i] = dataFile ( ( section + "/LinearParameters/Celerity1" ).data(), 1. );
212  M_celerity2[i] = dataFile ( ( section + "/LinearParameters/Celerity2" ).data(), 1. );
213  M_celerity1LeftEigenvector1[i] = dataFile ( ( section + "/LinearParameters/LeftEigenvector11" ).data(), 1. );
214  M_celerity1LeftEigenvector2[i] = dataFile ( ( section + "/LinearParameters/LeftEigenvector12" ).data(), 0. );
215  M_celerity2LeftEigenvector1[i] = dataFile ( ( section + "/LinearParameters/LeftEigenvector21" ).data(), 0. );
216  M_celerity2LeftEigenvector2[i] = dataFile ( ( section + "/LinearParameters/LeftEigenvector22" ).data(), 1. );
217  M_source10[i] = dataFile ( ( section + "/LinearParameters/Source10" ).data(), 0. );
218  M_source20[i] = dataFile ( ( section + "/LinearParameters/Source20" ).data(), 0. );
219  M_source11[i] = dataFile ( ( section + "/LinearParameters/Source11" ).data(), 0. );
220  M_source12[i] = dataFile ( ( section + "/LinearParameters/Source12" ).data(), 0. );
221  M_source21[i] = dataFile ( ( section + "/LinearParameters/Source21" ).data(), 0. );
222  M_source22[i] = dataFile ( ( section + "/LinearParameters/Source22" ).data(), 0. );
223  }
224 
225  break;
226 
227  case linear:
228 
229  linearInterpolation ( M_thickness, dataFile, section + "/PhysicalParameters/thickness", 0. );
230  linearInterpolation ( M_viscoelasticCoefficient, dataFile, section + "/PhysicalWall/ViscoelasticCoefficient", 0. );
231 
232  linearInterpolation ( M_area0, dataFile, section + "/PhysicalParameters/Area0", M_PI );
233  linearInterpolation ( M_alpha, dataFile, section + "/PhysicalParameters/AlphaCoriolis", 1. / M_robertsonCorrection, true );
234  linearInterpolation ( M_beta0, dataFile, section + "/PhysicalParameters/Beta0", 1.e6 );
235  linearInterpolation ( M_beta1, dataFile, section + "/PhysicalParameters/Beta1", 0.5 );
236 
237  linearInterpolation ( M_flux11, dataFile, section + "/PhysicalParameters/Flux11", 1. );
238  linearInterpolation ( M_flux12, dataFile, section + "/PhysicalParameters/Flux12", 0. );
239  linearInterpolation ( M_flux21, dataFile, section + "/PhysicalParameters/Flux21", 0. );
240  linearInterpolation ( M_flux22, dataFile, section + "/PhysicalParameters/Flux22", 1. );
241  linearInterpolation ( M_celerity1, dataFile, section + "/PhysicalParameters/Celerity1", 1. );
242  linearInterpolation ( M_celerity2, dataFile, section + "/PhysicalParameters/Celerity2", 1. );
243  linearInterpolation ( M_celerity1LeftEigenvector1, dataFile, section + "/PhysicalParameters/LeftEigenvector11", 1. );
244  linearInterpolation ( M_celerity1LeftEigenvector2, dataFile, section + "/PhysicalParameters/LeftEigenvector12", 0. );
245  linearInterpolation ( M_celerity2LeftEigenvector1, dataFile, section + "/PhysicalParameters/LeftEigenvector21", 0. );
246  linearInterpolation ( M_celerity2LeftEigenvector2, dataFile, section + "/PhysicalParameters/LeftEigenvector22", 1. );
247  linearInterpolation ( M_source10, dataFile, section + "/PhysicalParameters/Source10", 0. );
248  linearInterpolation ( M_source20, dataFile, section + "/PhysicalParameters/Source20", 0. );
249  linearInterpolation ( M_source11, dataFile, section + "/PhysicalParameters/Source11", 0. );
250  linearInterpolation ( M_source12, dataFile, section + "/PhysicalParameters/Source12", 0. );
251  linearInterpolation ( M_source21, dataFile, section + "/PhysicalParameters/Source21", 0. );
252  linearInterpolation ( M_source22, dataFile, section + "/PhysicalParameters/Source22", 0. );
253 
254  break;
255 
256  case pointwise:
257 
258  for ( UInt i = 0; i < M_meshPtr->numPoints() ; ++i )
259  {
260  // Physical Parameters
261  M_thickness[i] = dataFile ( ( section + "/PhysicalParameters/thickness" ).data(), 0., i );
262  M_viscoelasticCoefficient[i] = dataFile ( ( section + "/PhysicalWall/ViscoelasticCoefficient" ).data(), 0., i );
263 
264  M_area0[i] = dataFile ( ( section + "/PhysicalParameters/Area0" ).data(), M_PI, i );
265  M_alpha[i] = dataFile ( ( section + "/PhysicalParameters/AlphaCoriolis" ).data(), 1. / M_robertsonCorrection, i );
266  M_beta0[i] = dataFile ( ( section + "/PhysicalParameters/Beta0" ).data(), 1.e6, i );
267  M_beta1[i] = dataFile ( ( section + "/PhysicalParameters/Beta1" ).data(), 0.5, i );
268 
269  // Linear Parameters
270  M_flux11[i] = dataFile ( ( section + "/LinearParameters/Flux11" ).data(), 1., i );
271  M_flux12[i] = dataFile ( ( section + "/LinearParameters/Flux12" ).data(), 0., i );
272  M_flux21[i] = dataFile ( ( section + "/LinearParameters/Flux21" ).data(), 0., i );
273  M_flux22[i] = dataFile ( ( section + "/LinearParameters/Flux22" ).data(), 1., i );
274  M_celerity1[i] = dataFile ( ( section + "/LinearParameters/Celerity1" ).data(), 1., i );
275  M_celerity2[i] = dataFile ( ( section + "/LinearParameters/Celerity2" ).data(), 1., i );
276  M_celerity1LeftEigenvector1[i] = dataFile ( ( section + "/LinearParameters/LeftEigenvector11" ).data(), 1., i );
277  M_celerity1LeftEigenvector2[i] = dataFile ( ( section + "/LinearParameters/LeftEigenvector12" ).data(), 0., i );
278  M_celerity2LeftEigenvector1[i] = dataFile ( ( section + "/LinearParameters/LeftEigenvector21" ).data(), 0., i );
279  M_celerity2LeftEigenvector2[i] = dataFile ( ( section + "/LinearParameters/LeftEigenvector22" ).data(), 1., i );
280  M_source10[i] = dataFile ( ( section + "/LinearParameters/Source10" ).data(), 0., i );
281  M_source20[i] = dataFile ( ( section + "/LinearParameters/Source20" ).data(), 0., i );
282  M_source11[i] = dataFile ( ( section + "/LinearParameters/Source11" ).data(), 0., i );
283  M_source12[i] = dataFile ( ( section + "/LinearParameters/Source12" ).data(), 0., i );
284  M_source21[i] = dataFile ( ( section + "/LinearParameters/Source21" ).data(), 0., i );
285  M_source22[i] = dataFile ( ( section + "/LinearParameters/Source22" ).data(), 0., i );
286  }
287 
288  break;
289 
290  default:
291 
292  std::cout << "Warning: distributionLaw \"" << distributionLaw << "\"not available!" << std::endl;
293 
294  break;
295  }
296 
297  updateCoefficients();
298 }
299 
300 void
301 OneDFSIData::oldStyleSetup ( const GetPot& dataFile, const std::string& section )
302 {
303  std::cerr << std::endl << "Warning: OneDFSIData::oldStyleSetup( ... ) is deprecated!" << std::endl
304  << " You should use OneDFSIData::setup( ... ) instead!" << std::endl;
305 
306  // Model Type
307  M_physicsType = OneDFSI::physicsMap[ dataFile ( ( section + "/Model/PhysicsType" ).data(), "OneD_1DLinearPhysics" ) ];
308  M_fluxType = OneDFSI::fluxMap[ dataFile ( ( section + "/Model/FluxType" ).data(), "OneD_1DLinearFlux" ) ];
309  M_sourceType = OneDFSI::sourceMap[ dataFile ( ( section + "/Model/SourceType" ).data(), "OneD_1DLinearSource" ) ];
310 
311  // If data time has not been set
312  if ( !M_timeDataPtr.get() )
313  {
314  M_timeDataPtr.reset ( new time_Type ( dataFile, (section + "/time_discretization" ).data() ) );
315  }
316 
317  // Mesh setup - Space Discretization
318  Real length = dataFile ( ( section + "/discretization/x_right" ).data(), 1. ) -
319  dataFile ( ( section + "/discretization/x_left" ).data(), 0. );
320 
321  Real numberOfElements = dataFile ( ( section + "/discretization/nb_elem" ).data(), 10 );
322  regularMesh1D ( *M_meshPtr, 1, numberOfElements, false, length, 0 );
323 
324  //std::cout << " 1D- Mesh nodes: " << M_meshPtr->numPoints() << std::endl;
325  //std::cout << " 1D- Mesh elements: " << M_meshPtr->numElements() << std::endl;
326 
327  // Physical Wall Model
328  M_inertialWall = dataFile ( ( section + "/miscellaneous/inertial_wall" ).data(), false );
329  M_viscoelasticWall = dataFile ( ( section + "/miscellaneous/viscoelastic_wall" ).data(), false );
330  M_longitudinalWall = dataFile ( ( section + "/miscellaneous/longitudinal_wall" ).data(), false );
331 
332  // Miscellaneous
333  M_postprocessingDirectory = dataFile ( ( section + "/miscellaneous/post_dir" ).data(), "./" );
334  M_postprocessingFile = dataFile ( ( section + "/miscellaneous/post_file" ).data(), "sol" );
335  M_inertialWall = dataFile ( ( section + "/miscellaneous/inertial_wall" ).data(), false );
336  M_viscoelasticWall = dataFile ( ( section + "/miscellaneous/viscoelastic_wall" ).data(), false );
337  M_CFLmax = dataFile ( ( section + "/miscellaneous/CFLmax" ).data(), std::sqrt (3.) / 3. );
338 
339  if ( M_CFLmax > std::sqrt (3.) / 3. )
340  {
341  std::cout << "!!! WARNING: CFLmax greater than the theoretical value (see MOX21, eq. 1.47) - CONVERGENCE NOT GUARANTEED !!!" << std::endl;
342  }
343 
344  // Jacobian perturbation
345  M_jacobianPerturbationArea = dataFile ( ( section + "JacobianPerturbation/deltaArea" ).data(), 0.001 );
346  M_jacobianPerturbationFlowRate = dataFile ( ( section + "JacobianPerturbation/deltaFlowRate" ).data(), 0.001 );
347  M_jacobianPerturbationStress = dataFile ( ( section + "JacobianPerturbation/deltaStress" ).data(), 1.0 );
348 
349  // Physical Parameters
350  M_computeCoefficients = dataFile ( ( section + "/parameters/use_physical_values" ).data(), false );
351  M_powerLawCoefficient = dataFile ( ( section + "/PhysicalParameters/PowerlawCoefficient" ).data(), 2 );
352 
353  M_density = dataFile ( ( section + "/1d_physics/density" ).data(), 1. );
354  M_viscosity = dataFile ( ( section + "/1d_physics/viscosity" ).data(), 0.035 );
355 
356  M_densityWall = dataFile ( ( section + "/PhysicalParameters/densityWall" ).data(), 1. );
357  M_thickVessel = dataFile ( ( section + "/PhysicalParameters/thickVessel" ).data(), false );
358  M_young = dataFile ( ( section + "/1d_physics/young" ).data(), 4.0E6 );
359  M_poisson = dataFile ( ( section + "/1d_physics/ksi" ).data(), 0.5 );
360 
361  M_externalPressure = dataFile ( ( section + "/PhysicalParameters/externalPressure" ).data(), 0. );
362  M_venousPressure = dataFile ( ( section + "/PhysicalParameters/venousPressure" ).data(), 0. );
363 
364  M_inertialModulus = dataFile ( ( section + "/PhysicalParameters/coeffA" ).data(), 0. );
365  M_robertsonCorrection = dataFile ( ( section + "/PhysicalParameters/RobertsonCorrection" ).data(), 1. );
366 
367  M_friction = dataFile ( ( section + "/parameters/Kr" ).data(), 1. );
368  M_computeCoefficients = dataFile ( ( section + "/parameters/use_physical_values" ).data(), false );
369 
370  // Reset all the containers
371  resetContainers();
372 
373  // Linear Parameters
374  for ( UInt i = 0; i < M_meshPtr->numPoints() ; ++i )
375  {
376  M_viscoelasticCoefficient[i] = dataFile ( ( section + "/PhysicalParameters/gamma" ).data(), 0. );
377 
378  // Physical Parameters
379  if ( M_computeCoefficients )
380  {
381  M_area0[i] = OneDFSI::pow20 ( dataFile ( ( section + "/1d_physics/radius" ).data(), 0.5 ), 2 ) * M_PI;
382  }
383  else
384  {
385  M_area0[i] = dataFile ( ( section + "/parameters/Area0" ).data(), M_PI );
386  }
387  M_beta0[i] = dataFile ( ( section + "/parameters/beta0" ).data(), 1.e6 );
388  M_beta1[i] = dataFile ( ( section + "/parameters/beta1" ).data(), 0.5 );
389  M_thickness[i] = dataFile ( ( section + "/1d_physics/thickness" ).data(), 0. );
390 
391  M_alpha[i] = dataFile ( ( section + "/parameters/alphaCor" ).data(), 1. / M_robertsonCorrection );
392 
393  // Linear Parameters
394  M_flux11[i] = dataFile ( ( section + "/LinearParameters/Flux11" ).data(), 1. );
395  M_flux12[i] = dataFile ( ( section + "/LinearParameters/Flux12" ).data(), 0. );
396  M_flux21[i] = dataFile ( ( section + "/LinearParameters/Flux21" ).data(), 0. );
397  M_flux22[i] = dataFile ( ( section + "/LinearParameters/Flux22" ).data(), 1. );
398  M_celerity1[i] = dataFile ( ( section + "/LinearParameters/Celerity1" ).data(), 1. );
399  M_celerity2[i] = dataFile ( ( section + "/LinearParameters/Celerity2" ).data(), 1. );
400  M_celerity1LeftEigenvector1[i] = dataFile ( ( section + "/LinearParameters/LeftEigenvector11" ).data(), 1. );
401  M_celerity1LeftEigenvector2[i] = dataFile ( ( section + "/LinearParameters/LeftEigenvector12" ).data(), 0. );
402  M_celerity2LeftEigenvector1[i] = dataFile ( ( section + "/LinearParameters/LeftEigenvector21" ).data(), 0. );
403  M_celerity2LeftEigenvector2[i] = dataFile ( ( section + "/LinearParameters/LeftEigenvector22" ).data(), 1. );
404  M_source10[i] = dataFile ( ( section + "/LinearParameters/Source10" ).data(), 0. );
405  M_source20[i] = dataFile ( ( section + "/LinearParameters/Source20" ).data(), 0. );
406  M_source11[i] = dataFile ( ( section + "/LinearParameters/Source11" ).data(), 0. );
407  M_source12[i] = dataFile ( ( section + "/LinearParameters/Source12" ).data(), 0. );
408  M_source21[i] = dataFile ( ( section + "/LinearParameters/Source21" ).data(), 0. );
409  M_source22[i] = dataFile ( ( section + "/LinearParameters/Source22" ).data(), 0. );
410  }
411 
412  updateCoefficients();
413 }
414 
415 void
416 OneDFSIData::updateCoefficients()
417 {
418  if ( M_computeCoefficients )
419  {
420  // PowerlawProfile: s(r) = (1+2/gamma)*(1-r^gamma)
421  Real radius (1.); //std::sqrt( M_Area0[i] / M_PI );
422 
423  Real profileIntegral = ( 1 + 2. / M_powerLawCoefficient ) * ( 1 + 2. / M_powerLawCoefficient ) *
424  ( radius * radius / 2 +
425  std::pow (radius, 2 * M_powerLawCoefficient + 2) / (2 * M_powerLawCoefficient + 2 ) -
426  2 * std::pow (radius, M_powerLawCoefficient + 2) / ( M_powerLawCoefficient + 2 )
427  );
428 
429  // Compute Friction Coefficient: Kr = -2*pi*mu/rho*s'(R)
430  M_friction = 2 * M_PI * M_viscosity / M_density * ( M_powerLawCoefficient + 2 ) * std::pow ( radius, M_powerLawCoefficient - 1 );
431 
432  for ( UInt i = 0; i < M_meshPtr->numPoints(); ++i )
433  {
434  // Compute Coriolis Coefficient: Alpha = 2*pi/Area0*Int(s(r)^2)
435  M_alpha[i] = 2 / ( radius * radius ) * profileIntegral;
436 
437  // Compute Beta0
438  if ( M_thickVessel ) // see Amadori, Ferrari, Formaggia (MOX report 86)
439  {
440  M_beta0[i] = - M_thickness[i] * M_young * std::sqrt (M_PI) / ( std::sqrt ( M_area0[i] ) * ( (1 - M_poisson * M_poisson )
441  + M_poisson * ( 1 + M_poisson ) * ( M_thickness[i] * std::sqrt (M_PI) / std::sqrt ( M_area0[i] ) ) ) );
442  M_beta1[i] = - 0.5;
443  }
444  else
445  {
446  M_beta0[i] = M_thickness[i] * std::sqrt ( M_PI / M_area0[i] ) * M_young / (1 - M_poisson * M_poisson);
447  M_beta1[i] = 0.5;
448  }
449 
450  // Compute Viscoelastic Coefficient: gamma = h * E / ( 1 - nu^2 ) * T * tan(phi) / ( 4 * \sqrt(pi))
451  M_viscoelasticCoefficient[i] = M_thickness[i] * M_young / (1 - M_poisson * M_poisson) * M_viscoelasticPeriod *
452  std::tan ( M_viscoelasticAngle ) / ( 4 * std::sqrt (M_PI) );
453  }
454  }
455 
456  // Compute the derivatives of the coefficients
457  computeSpatialDerivatives();
458 }
459 
460 //void
461 //OneDFSIData::initializeLinearParameters()
462 //{
463 // for ( UInt indz=0; indz < M_meshPtr->numPoints(); ++indz )
464 // {
465 // M_celerity1[indz] = std::sqrt( M_beta0[indz] * M_beta1[indz] / M_density );
466 // M_flux21[indz] = M_celerity1[indz] * M_celerity1[indz];
467 // M_source22[indz] = M_friction / M_area0(indz);
468 // }
469 //
470 // M_celerity2 = - M_celerity1;
471 //
472 // M_celerity1LeftEigenvector1 = M_celerity1;
473 // M_celerity1LeftEigenvector2 = scalarVector_Type(M_meshPtr->numPoints(), 1.);
474 // M_celerity2LeftEigenvector1 = M_celerity2;
475 // M_celerity2LeftEigenvector2 = scalarVector_Type(M_meshPtr->numPoints(), 1.);
476 //
477 // M_flux11 = ZeroVector(M_meshPtr->numPoints());
478 // M_flux12 = scalarVector_Type(M_meshPtr->numPoints(), 1.);
479 // M_flux22 = ZeroVector(M_meshPtr->numPoints());
480 //
481 // M_source10 = ZeroVector(M_meshPtr->numPoints());
482 // M_source11 = ZeroVector(M_meshPtr->numPoints());
483 // M_source12 = ZeroVector(M_meshPtr->numPoints());
484 //
485 // M_source20 = ZeroVector(M_meshPtr->numPoints());
486 // M_source21 = ZeroVector(M_meshPtr->numPoints());
487 //}
488 
489 // TODO Many changes to be done in these two methods:
490 // 1) variables name is very bad
491 // 2) do loop with explicit i++ is bad, use for loop instead
492 // 3) rename indices following other classes (iNode, iElement, etc..)
493 // 4) probably 3/4 of the code can be shared between the left and right method instead of duplications
494 //void
495 //OneDFSIData::stiffenVesselLeft( const Real& xl, const Real& xr,
496 // const Real& factor, const Real& alpha,
497 // const Real& delta, const Real& n,
498 // const Real& minDeltaX, const UInt& yesAdaptive )
499 //{
500 // if ( yesAdaptive )
501 // {
502 // UInt iz=0, alpha_iz;
503 // Real ratio, n_elem_delta,n_elem_r,n_elem_l;
504 //
505 // alpha_iz = static_cast<Int>( std::floor( (alpha - delta / 2 ) / minDeltaX + 0.5 ) ) + ( ( numberOfElements() - 1 ) -
506 // static_cast<Int>( std::floor( ( xr - ( alpha + delta / 2 ) ) / minDeltaX + 0.5 ) ) -
507 // static_cast<Int>( std::floor( ( alpha - delta / 2 ) / minDeltaX + 0.5 ) ) ) / 2;
508 //
509 // n_elem_r = ( ( numberOfElements() - 1 ) - alpha_iz ) -
510 // static_cast<Int>( std::floor( ( xr - ( alpha + delta / 2 ) ) / minDeltaX + 0.5 ) );
511 //
512 // n_elem_l = alpha_iz -
513 // static_cast<Int>( std::floor( ( alpha - delta / 2 ) / minDeltaX + 0.5 ) );
514 //
515 // n_elem_delta = static_cast<Real>( numberOfElements() - 1 ) / ( xr - xl ) * delta;
516 //
517 // Real x_current,deltax,deltax_adaptive,deltax_uniform;
518 // x_current = alpha;
519 // do
520 // {
521 // ratio = ( alpha + delta / 2 - x_current ) / delta;
522 //
523 // M_dBeta0dz[alpha_iz + iz] = M_beta0[alpha_iz + iz] * ( factor * (- n / delta) * ( std::pow(2,(n-1)) * std::pow(ratio, (n-1)) ) );
524 // M_dBeta0dz[alpha_iz - iz] = M_dBeta0dz[alpha_iz + iz];
525 //
526 // M_beta0[alpha_iz + iz] = M_beta0[alpha_iz + iz] * ( 1 + factor * ( std::pow(2,(n-1)) * std::pow(ratio,n) ) );
527 // M_beta0[alpha_iz - iz] = M_beta0[alpha_iz + iz] / ( 1 + factor * ( std::pow(2,(n-1)) * std::pow(ratio,n) ) ) * ( 1 + factor * ( 1 - ( std::pow(2,(n-1)) * std::pow(ratio,n) ) ) );
528 //
529 // deltax_adaptive = ( -1 / n_elem_delta ) * ( 1 / ( ( -n / delta) * std::pow(2,(n-1) ) * std::pow( ratio , (n-1) ) ) );
530 // deltax_uniform = ( ( alpha + delta / 2) - x_current ) / ( n_elem_r - iz );
531 //
532 // ++iz;
533 // deltax = ( ( deltax_adaptive < deltax_uniform ) && ( iz < n_elem_r ) ) ? deltax_adaptive : deltax_uniform;
534 // ASSERT_PRE( deltax > 0 , "The left point is on the right..." );
535 // x_current += deltax;
536 // }
537 // while ( ( x_current < ( alpha + delta/2 ) ) && ( ( alpha_iz - ( iz - 1 ) ) > 0) );
538 //
539 // if ( ( alpha_iz - ( iz - 1 ) ) > 0 )
540 // {
541 // do
542 // {
543 // M_beta0[alpha_iz - iz] = M_beta0[alpha_iz - iz] * ( 1 + factor );
544 // ++iz;
545 // }
546 // while ( (alpha_iz - ( iz - 1 ) ) > 0 );
547 // }
548 // else
549 // std::cout << "[stiffenVesselRight] error! out of left boundary" << std::endl;
550 // }
551 //
552 // else
553 // {
554 // UInt iz=0;
555 // Real ratio, x_current=xl, deltax;
556 //
557 // deltax = ( xr - xl ) / static_cast<Real>(numberOfElements() - 1);
558 //
559 // while ( ( x_current < ( alpha - delta / 2 ) ) && ( iz < numberOfElements() ) )
560 // {
561 // M_beta0[iz] = M_beta0[iz] * ( 1 + factor );
562 //
563 // ++iz;
564 // x_current+=deltax;
565 // }
566 //
567 // while ( (x_current < alpha) && (iz < numberOfElements()) )
568 // {
569 // ratio = ( x_current - alpha - delta / 2 ) / delta;
570 //
571 // M_dBeta0dz[iz] = M_beta0[iz] * ( factor * (- n / delta ) * ( std::pow(2,(n-1)) * std::pow(ratio,(n-1) ) ) );
572 // M_beta0[iz] = M_beta0[iz] * ( 1 + factor * ( 1 - std::pow(2,(n-1)) * std::pow(ratio,n) ) );
573 //
574 // ++iz;
575 // x_current+=deltax;
576 // }
577 //
578 // while ( ( x_current < ( alpha + delta / 2 ) ) && (iz < numberOfElements()) )
579 // {
580 // ratio = ( alpha + delta / 2 - x_current ) / delta;
581 //
582 // M_dBeta0dz[iz] = M_beta0[iz] * ( factor * ( -n / delta) * ( std::pow(2,(n-1)) * std::pow(ratio,(n-1) ) ) );
583 // M_beta0[iz] = M_beta0[iz] * ( 1 + factor * ( std::pow(2,(n-1)) * std::pow(ratio,n) ) );
584 //
585 // ++iz;
586 // x_current += deltax;
587 // }
588 // }
589 //}
590 
591 //void
592 //OneDFSIData::stiffenVesselRight( const Real& xl, const Real& xr,
593 // const Real& factor, const Real& alpha,
594 // const Real& delta, const Real& n,
595 // const Real& minDeltaX, const UInt& yesAdaptive )
596 //{
597 // if ( yesAdaptive )
598 // {
599 // Real ratio, n_elem_delta,n_elem_r,n_elem_l;
600 //
601 // UInt iz=0, alpha_iz;
602 //
603 // alpha_iz = static_cast<Int>( std::floor( ( alpha - delta / 2 ) / minDeltaX + 0.5 ) ) + ( (numberOfElements() - 1) -
604 // static_cast<Int>( std::floor( ( xr - ( alpha + delta / 2 ) ) / minDeltaX + 0.5 ) ) -
605 // static_cast<Int>( std::floor( ( alpha - delta / 2 ) / minDeltaX + 0.5 ) ) ) / 2;
606 //
607 // n_elem_delta = static_cast<Real>( numberOfElements() - 1 ) / ( xr - xl ) * delta;
608 //
609 // n_elem_r = ( ( numberOfElements() - 1 ) - alpha_iz ) -
610 // static_cast<Int>( std::floor( ( xr - ( alpha + delta / 2 ) ) / minDeltaX + 0.5 ) );
611 //
612 // n_elem_l = alpha_iz -
613 // static_cast<Int>( std::floor( ( alpha - delta / 2 ) / minDeltaX + 0.5 ) );
614 //
615 // Real x_current,deltax,deltax_adaptive,deltax_uniform;
616 // x_current = alpha;
617 // do
618 // {
619 // ratio = ( alpha + delta / 2 - x_current ) / delta;
620 //
621 // M_dBeta0dz[alpha_iz + iz] = M_beta0[alpha_iz + iz] * ( factor * ( n / delta) * ( std::pow(2,(n-1)) * std::pow(ratio,(n-1)) ) );
622 // M_dBeta0dz[alpha_iz - iz] = M_dBeta0dz[alpha_iz + iz];
623 //
624 // M_beta0[alpha_iz + iz] = M_beta0[alpha_iz + iz] * ( 1 + factor * ( 1 - ( std::pow(2,(n-1)) * std::pow(ratio,n) ) ) );
625 // M_beta0[alpha_iz - iz] = M_beta0[alpha_iz + iz] / ( 1 + factor * ( 1 - ( std::pow(2,(n-1)) * std::pow(ratio,n) ) ) ) * ( 1 + factor * ( std::pow(2,(n-1)) * std::pow(ratio,n) ) );
626 //
627 // deltax_adaptive = ( -1 / n_elem_delta ) * ( 1 / ( ( -n / delta) * std::pow(2,(n-1) ) * std::pow( ratio , (n-1) ) ) );
628 // deltax_uniform = ( ( alpha + delta / 2) - x_current ) / ( n_elem_r - iz );
629 //
630 // ++iz;
631 // deltax = ( ( deltax_adaptive < deltax_uniform ) && ( iz < n_elem_r) ) ? deltax_adaptive : deltax_uniform;
632 // ASSERT_PRE( deltax > 0 , "The left point is on the right..." );
633 // x_current += deltax;
634 //
635 // }
636 // while ( x_current < ( alpha + delta / 2 ) && ( ( alpha_iz - ( iz - 1 ) ) > 0) );
637 //
638 // if ( ( alpha_iz + iz ) <= ( numberOfElements() -1 ) )
639 // {
640 // do
641 // {
642 // M_beta0[alpha_iz + iz] = M_beta0[alpha_iz + iz] * ( 1 + factor );
643 // ++iz;
644 // }
645 // while ( ( alpha_iz + iz - 1 ) < ( numberOfElements() -1 ) );
646 //
647 // }
648 // else
649 // std::cout << "\n[stiffenVesselRight] error! out of right boundary" << std::endl;
650 // }
651 // else
652 // {
653 // UInt iz = numberOfElements()-1;
654 // Real ratio, x_current=xr, deltax;
655 //
656 // deltax = ( xr - xl ) / static_cast<Real>(numberOfElements() - 1 );
657 //
658 // while ( ( x_current > ( alpha + delta / 2 ) ) && ( ( iz + 1 ) > 0 ) )
659 // {
660 // M_beta0[iz] = M_beta0[iz] * ( 1 + factor );
661 //
662 // --iz;
663 // x_current -= deltax;
664 // }
665 //
666 // while ( ( x_current > alpha ) && ( ( iz + 1 ) > 0 ) )
667 // {
668 // ratio=( ( ( alpha + delta / 2 ) - x_current ) / delta );
669 //
670 // M_dBeta0dz[iz] = M_beta0[iz] * ( factor * ( n / delta) * ( std::pow(2,(n-1)) * std::pow(ratio,(n-1)) ) );
671 // M_beta0[iz] = M_beta0[iz] * ( 1 + factor * ( 1 - std::pow(2,(n-1)) * std::pow(ratio,n) ) );
672 //
673 // --iz;
674 // x_current -= deltax;
675 // }
676 //
677 // while ( ( x_current > ( alpha - delta / 2 ) ) && ( ( iz + 1 ) > 0 ) )
678 // {
679 // ratio = ( x_current - alpha - delta / 2 ) / delta;
680 //
681 // M_dBeta0dz[iz] = M_beta0[iz] * ( factor * ( n / delta) * ( std::pow(2,(n-1)) * std::pow(ratio,(n-1)) ) );
682 // M_beta0[iz] = M_beta0[iz] * ( 1 + factor * ( std::pow(2,(n-1)) * std::pow(ratio,n) ) );
683 //
684 // --iz;
685 // x_current -= deltax;
686 // }
687 // }
688 //}
689 
690 void
691 OneDFSIData::showMe ( std::ostream& output ) const
692 {
693  //output << std::scientific << std::setprecision(15);
694 
695  // Model
696  output << "\n*** Values for data [Model]" << std::endl << std::endl;
697  output << "Physics Type = " << enum2String ( M_physicsType, OneDFSI::physicsMap ) << std::endl;
698  output << "Flux Type = " << enum2String ( M_fluxType, OneDFSI::fluxMap ) << std::endl;
699  output << "Source Type = " << enum2String ( M_sourceType, OneDFSI::sourceMap ) << std::endl;
700 
701  // Time
702  output << "\n*** Values for data [time_discretization]" << std::endl << std::endl;
703  M_timeDataPtr->showMe ( output );
704 
705  // Space Discretization
706  output << "\n*** Values for data [space_discretization]" << std::endl << std::endl;
707  output << "Length = " << length() << std::endl;
708  output << "NumberOfElements = " << M_meshPtr->numElements() << std::endl;
709 
710  // Physical Wall Model
711  output << "\n*** Values for data [PhysicalWallModel]" << std::endl << std::endl;
712  output << "Viscoelastic wall = " << M_viscoelasticWall << std::endl;
713  output << "Viscoelastic angle = " << M_viscoelasticAngle << std::endl;
714  output << "Viscoelastic period = " << M_viscoelasticPeriod << std::endl;
715  output << "Viscoelastic coeff. = " << M_viscoelasticCoefficient << std::endl;
716  output << "Inertial wall = " << M_inertialWall << std::endl;
717  output << "Wall density = " << M_densityWall << std::endl;
718  output << "Inertial modulus = " << M_inertialModulus << std::endl;
719  output << "Longitudinal wall = " << M_longitudinalWall << std::endl;
720 
721  // Miscellaneous
722  output << "\n*** Values for data [miscellaneous]" << std::endl << std::endl;
723  output << "Postprocessing Dir. = " << M_postprocessingDirectory << std::endl;
724  output << "Postprocessing File = " << M_postprocessingFile << std::endl;
725  output << "Maximum admissible CFL = " << M_CFLmax << std::endl;
726 
727  // Jacobian perturbation
728  output << "Jacobian perturbation Area = " << M_jacobianPerturbationArea << std::endl;
729  output << "Jacobian perturbation Flow Rate = " << M_jacobianPerturbationFlowRate << std::endl;
730  output << "Jacobian perturbation Stress = " << M_jacobianPerturbationStress << std::endl;
731 
732  // Physical Parameters
733  output << "\n*** Values for data [PhysicalParameters]" << std::endl << std::endl;
734  output << "Compute Coefficients = " << M_computeCoefficients << std::endl;
735  output << "Powerlaw Coefficient = " << M_powerLawCoefficient << std::endl;
736  output << "Fluid density = " << M_density << std::endl;
737  output << "Fluid dyn. viscosity = " << M_viscosity << std::endl;
738  output << "Thick vessel = " << M_thickVessel << std::endl;
739  output << "Young modulus = " << M_young << std::endl;
740  output << "Poisson = " << M_poisson << std::endl;
741  output << "External pressure = " << M_externalPressure << std::endl;
742  output << "Venous pressure = " << M_venousPressure << std::endl;
743  output << "Robertson correction = " << M_robertsonCorrection << std::endl;
744 
745  output << "Area0 = " << M_area0 << std::endl;
746  output << "dArea0dz = " << M_dArea0dz << std::endl;
747  output << "Beta0 = " << M_beta0 << std::endl;
748  output << "dBeta0dz = " << M_dBeta0dz << std::endl;
749  output << "Beta1 = " << M_beta1 << std::endl;
750  output << "dBeta1dz = " << M_dBeta1dz << std::endl;
751 
752  output << "Alpha (Coriolis) = " << M_alpha << std::endl;
753  output << "dAlpha (Coriolis) = " << M_dAlphadz << std::endl;
754 
755  output << "Thickness = " << M_thickness << std::endl;
756  output << "Friction = " << M_friction << std::endl;
757 
758  // Linear Parameters
759  output << "\n*** Values for data [LinearParameters]" << std::endl << std::endl;
760  output << "Flux11 = " << M_flux11 << std::endl;
761  output << "Flux12 = " << M_flux12 << std::endl;
762  output << "Flux21 = " << M_flux21 << std::endl;
763  output << "Flux22 = " << M_flux22 << std::endl;
764  output << "Celerity1 = " << M_celerity1 << std::endl;
765  output << "Celerity2 = " << M_celerity2 << std::endl;
766  output << "Eigenvector11 = " << M_celerity1LeftEigenvector1 << std::endl;
767  output << "Eigenvector12 = " << M_celerity1LeftEigenvector2 << std::endl;
768  output << "Eigenvector21 = " << M_celerity2LeftEigenvector1 << std::endl;
769  output << "Eigenvector22 = " << M_celerity2LeftEigenvector2 << std::endl;
770  output << "Source10 = " << M_source10 << std::endl;
771  output << "Source20 = " << M_source20 << std::endl;
772  output << "Source11 = " << M_source11 << std::endl;
773  output << "Source12 = " << M_source12 << std::endl;
774  output << "Source21 = " << M_source21 << std::endl;
775  output << "Source22 = " << M_source22 << std::endl;
776 }
777 
778 // ===================================================
779 // Get Methods - Physical Parameters
780 // ===================================================
781 const Real&
782 OneDFSIData::robertsonCorrection() const
783 {
784  if ( M_robertsonCorrection != 1. )
785  {
786  std::cout << "!!! WARNING: Robertson corretion has not been checked in this version of the code !!!" << std::endl;
787  }
788 
789  return M_robertsonCorrection;
790 }
791 // ===================================================
792 // Private methods
793 // ===================================================
794 void
795 OneDFSIData::linearInterpolation ( scalarVector_Type& vector,
796  const GetPot& dataFile,
797  const std::string& quantity,
798  const Real& defaultValue,
799  const bool& isArea )
800 {
801  Real a = dataFile ( quantity.data(), defaultValue, 0 );
802  Real b = dataFile ( quantity.data(), a, 1 );
803 
804  Real xa = M_meshPtr->firstPoint().x();
805  Real xb = M_meshPtr->lastPoint().x();
806 
807  for ( UInt i (0) ; i < M_meshPtr->numPoints() ; ++i )
808  if ( isArea )
809  {
810  vector[i] = std::sqrt (a / M_PI) + ( std::sqrt (b / M_PI) - std::sqrt (a / M_PI) ) /
811  ( xb - xa ) * ( M_meshPtr->point ( i ).x() - xa );
812  vector[i] *= vector[i] * M_PI;
813  }
814  else
815  {
816  vector[i] = a + (b - a) / ( xb - xa ) * ( M_meshPtr->point ( i ).x() - xa );
817  }
818 
819  // linearInterpolation to disable tapering (when needed)
820  // for ( UInt i(0); i < M_meshPtr->numPoints() ; ++i )
821  // if ( isArea )
822  // vector[i] = (a + b + 2 * std::sqrt(a*b)) / 4;
823  // else
824  // vector[i] = (a + b) / 2;
825 }
826 
827 void
828 OneDFSIData::computeSpatialDerivatives()
829 {
830  for ( UInt iNode ( 0 ) ; iNode < M_meshPtr->numPoints() ; ++iNode )
831  {
832  M_dArea0dz[iNode] = computeSpatialDerivativeAtNode ( M_area0, iNode );
833  M_dBeta0dz[iNode] = computeSpatialDerivativeAtNode ( M_beta0, iNode );
834  M_dBeta1dz[iNode] = computeSpatialDerivativeAtNode ( M_beta1, iNode );
835  M_dAlphadz[iNode] = computeSpatialDerivativeAtNode ( M_alpha, iNode );
836  }
837 }
838 
839 void
840 OneDFSIData::resetContainers()
841 {
842  M_viscoelasticCoefficient.resize ( M_meshPtr->numPoints() );
843  M_thickness.resize ( M_meshPtr->numPoints() );
844 
845  M_area0.resize ( M_meshPtr->numPoints() );
846  M_beta0.resize ( M_meshPtr->numPoints() );
847  M_beta1.resize ( M_meshPtr->numPoints() );
848  M_alpha.resize ( M_meshPtr->numPoints() );
849 
850  M_dArea0dz.resize ( M_meshPtr->numPoints() );
851  M_dBeta0dz.resize ( M_meshPtr->numPoints() );
852  M_dBeta1dz.resize ( M_meshPtr->numPoints() );
853  M_dAlphadz.resize ( M_meshPtr->numPoints() );
854 
855  // Linear Parameters defined along the 1D model
856  M_flux11.resize ( M_meshPtr->numPoints() );
857  M_flux12.resize ( M_meshPtr->numPoints() );
858  M_flux21.resize ( M_meshPtr->numPoints() );
859  M_flux22.resize ( M_meshPtr->numPoints() );
860  M_celerity1.resize ( M_meshPtr->numPoints() );
861  M_celerity2.resize ( M_meshPtr->numPoints() );
862  M_celerity1LeftEigenvector1.resize ( M_meshPtr->numPoints() );
863  M_celerity1LeftEigenvector2.resize ( M_meshPtr->numPoints() );
864  M_celerity2LeftEigenvector1.resize ( M_meshPtr->numPoints() );
865  M_celerity2LeftEigenvector2.resize ( M_meshPtr->numPoints() );
866  M_source10.resize ( M_meshPtr->numPoints() );
867  M_source20.resize ( M_meshPtr->numPoints() );
868  M_source11.resize ( M_meshPtr->numPoints() );
869  M_source12.resize ( M_meshPtr->numPoints() );
870  M_source21.resize ( M_meshPtr->numPoints() );
871  M_source22.resize ( M_meshPtr->numPoints() );
872 }
873 
874 } // LifeV namespace
LifeV::OneDFSIData::M_postprocessingDirectory
std::string M_postprocessingDirectory
Miscellaneous.
Definition: OneDFSIData.hpp:953
LifeV::OneDFSIData::M_flux12
scalarVector_Type M_flux12
Definition: OneDFSIData.hpp:993
LifeV::OneDFSIData::M_flux21
scalarVector_Type M_flux21
Definition: OneDFSIData.hpp:994
LifeV::OneDFSIData::OneDFSIData
OneDFSIData()
Empty constructor.
Definition: OneDFSIData.cpp:51
LifeV::OneDFSIData::M_celerity1LeftEigenvector2
scalarVector_Type M_celerity1LeftEigenvector2
Definition: OneDFSIData.hpp:1003
LifeV::OneDFSIData::M_dAlphadz
scalarVector_Type M_dAlphadz
Definition: OneDFSIData.hpp:987
LifeV::OneDFSIData::M_alpha
scalarVector_Type M_alpha
Definition: OneDFSIData.hpp:981
LifeV::OneDFSIData::M_thickness
scalarVector_Type M_thickness
Definition: OneDFSIData.hpp:977
LifeV::OneDFSIData
OneDFSIData - Class which read and holds all the data for the One Dimensional Model Solver...
Definition: OneDFSIData.hpp:106
LifeV::OneDFSIData::M_dArea0dz
scalarVector_Type M_dArea0dz
Derivatives of main coefficients.
Definition: OneDFSIData.hpp:986
LifeV::OneDFSIData::M_celerity2LeftEigenvector2
scalarVector_Type M_celerity2LeftEigenvector2
Definition: OneDFSIData.hpp:1005
LifeV::OneDFSIData::M_flux11
scalarVector_Type M_flux11
Flux matrix.
Definition: OneDFSIData.hpp:992
LifeV::OneDFSIData::M_source21
scalarVector_Type M_source21
Definition: OneDFSIData.hpp:1012
LifeV::OneDFSIData::M_computeCoefficients
bool M_computeCoefficients
Physical Parameters.
Definition: OneDFSIData.hpp:963
LifeV::OneDFSIData::M_source12
scalarVector_Type M_source12
Definition: OneDFSIData.hpp:1011
LifeV::OneDFSIData::updateCoefficients
void updateCoefficients()
Update all the physical coefficients.
Definition: OneDFSIData.cpp:416
M_PI
#define M_PI
Definition: winmath.h:20
LifeV::OneDFSIData::M_source20
scalarVector_Type M_source20
Definition: OneDFSIData.hpp:1009
LifeV::OneDFSIData::M_jacobianPerturbationStress
Real M_jacobianPerturbationStress
Definition: OneDFSIData.hpp:960
ETCurrentFE::updateInverseJacobian
void updateInverseJacobian(const UInt &iQuadPt)
Definition: ETCurrentFE.cpp:405
LifeV::OneDFSIData::resetContainers
void resetContainers()
Reset all the containers.
Definition: OneDFSIData.cpp:840
LifeV::OneDFSIData::M_source11
scalarVector_Type M_source11
Definition: OneDFSIData.hpp:1010
LifeV::OneDFSIData::M_sourceType
OneDFSI::sourceTerm_Type M_sourceType
Definition: OneDFSIData.hpp:935
LifeV::OneDFSIData::M_celerity1LeftEigenvector1
scalarVector_Type M_celerity1LeftEigenvector1
Eigenvector for first and second eigenvalue.
Definition: OneDFSIData.hpp:1002
LifeV::OneDFSIData::M_viscoelasticCoefficient
scalarVector_Type M_viscoelasticCoefficient
Definition: OneDFSIData.hpp:945
LifeV::OneDFSIData::robertsonCorrection
const Real & robertsonCorrection() const
Get the Robertson correction coefficient (Not tested: maybe wrong in the code)
Definition: OneDFSIData.cpp:782
LifeV::OneDFSIData::M_area0
scalarVector_Type M_area0
Definition: OneDFSIData.hpp:980
LifeV::OneDFSIData::M_celerity2LeftEigenvector1
scalarVector_Type M_celerity2LeftEigenvector1
Definition: OneDFSIData.hpp:1004
LifeV::OneDFSIData::showMe
void showMe(std::ostream &output=std::cout) const
Initialize linear parameters (NOT WORKING)
Definition: OneDFSIData.cpp:691
LifeV::OneDFSIData::M_postprocessingFile
std::string M_postprocessingFile
full directory name (including path)
Definition: OneDFSIData.hpp:954
LifeV::OneDFSIData::M_beta0
scalarVector_Type M_beta0
Definition: OneDFSIData.hpp:982
LifeV::OneDFSIData::M_beta1
scalarVector_Type M_beta1
Definition: OneDFSIData.hpp:983
LifeV::OneDFSIData::computeSpatialDerivatives
void computeSpatialDerivatives()
Compute the derivatives of alpha, area0, beta0, and beta1 using centered differences.
Definition: OneDFSIData.cpp:828
LifeV::OneDFSIData::scalarVector_Type
ublas::vector< Real > scalarVector_Type
Definition: OneDFSIData.hpp:120
LifeV::OneDFSIData::M_flux22
scalarVector_Type M_flux22
Definition: OneDFSIData.hpp:995
LifeV::OneDFSIData::M_dBeta1dz
scalarVector_Type M_dBeta1dz
Definition: OneDFSIData.hpp:989
LifeV::OneDFSIData::M_celerity1
scalarVector_Type M_celerity1
Celerities of the linear problem (eigenvalues of the flux matrix)
Definition: OneDFSIData.hpp:998
LifeV::Real
double Real
Generic real data.
Definition: LifeV.hpp:175
LifeV::OneDFSIData::M_jacobianPerturbationArea
Real M_jacobianPerturbationArea
Jacobian perturbation.
Definition: OneDFSIData.hpp:958
LifeV::OneDFSIData::M_source22
scalarVector_Type M_source22
Definition: OneDFSIData.hpp:1013
LifeV::OneDFSIData::M_physicsType
OneDFSI::physicsType_Type M_physicsType
Model.
Definition: OneDFSIData.hpp:933
LifeV::OneDFSIData::M_timeDataPtr
timePtr_Type M_timeDataPtr
Data containers for time and mesh.
Definition: OneDFSIData.hpp:938
LifeV::OneDFSIData::M_dBeta0dz
scalarVector_Type M_dBeta0dz
Definition: OneDFSIData.hpp:988
LifeV::OneDFSIData::M_source10
scalarVector_Type M_source10
Source matrix.
Definition: OneDFSIData.hpp:1008
LifeV::OneDFSIData::M_jacobianPerturbationFlowRate
Real M_jacobianPerturbationFlowRate
Definition: OneDFSIData.hpp:959
LifeV::OneDFSIData::M_fluxType
OneDFSI::fluxTerm_Type M_fluxType
Definition: OneDFSIData.hpp:934
LifeV::OneDFSIData::M_CFLmax
Real M_CFLmax
output file name
Definition: OneDFSIData.hpp:955
LifeV::OneDFSIData::M_viscoelasticWall
bool M_viscoelasticWall
Physical Wall Model.
Definition: OneDFSIData.hpp:942
LifeV::OneDFSIData::M_celerity2
scalarVector_Type M_celerity2
Definition: OneDFSIData.hpp:999