Commit 8092e131 authored by Simone Deparis's avatar Simone Deparis Committed by Simone Rossi
Browse files

Astyle on EM branch

parent d660e7e0
......@@ -860,7 +860,7 @@ public:
@param ionicSize number of equation in the ionic model
*/
virtual void setup (std::string meshName, std::string meshPath, GetPot& dataFile,
short int ionicSize);
short int ionicSize);
//! create mass matrix
/*!
......@@ -1741,7 +1741,7 @@ void ElectroETAMonodomainSolver<Mesh, IonicModel>::setupStiffnessMatrix (
{
using namespace ExpressionAssembly;
auto I = value(M_identity);
auto I = value (M_identity);
auto f0 = value (spaceVectorial, *M_fiberPtr);
auto D = value (sigmat) * I + (value (sigmal) - value (sigmat) ) * outerProduct (f0, f0);
......
......@@ -168,7 +168,7 @@ public:
/*!
* @param parameterList Teuchos parameter list
*/
virtual void setup( Teuchos::ParameterList& parameterList) { }
virtual void setup ( Teuchos::ParameterList& parameterList) { }
//! returns the number of equations of the ionic model
/*!
......
......@@ -63,7 +63,7 @@ IonicAlievPanfilov::IonicAlievPanfilov() :
IonicAlievPanfilov::IonicAlievPanfilov ( Teuchos::ParameterList& parameterList ) :
super ( 2 )
{
setup ( parameterList );
setup ( parameterList );
M_restingConditions.at (0) = 0.0;
M_restingConditions.at (1) = 0.16;
......
......@@ -61,7 +61,7 @@ IonicFitzHughNagumo::IonicFitzHughNagumo() :
IonicFitzHughNagumo::IonicFitzHughNagumo ( Teuchos::ParameterList& parameterList ) :
super ( 2 )
{
setup ( parameterList );
setup ( parameterList );
M_restingConditions.at (0) = 1e-8;
M_restingConditions.at (1) = 0.3;
......
......@@ -781,26 +781,26 @@ int main (int argc, char** argv)
boost::shared_ptr<HeavisideFct> H (new HeavisideFct);
BOOST_AUTO_TPL (deformationGradientTensor, ( grad( electrodETFESpace, *emDisp, 0) + value( solid.activeMaterial()-> identity() ) ));
BOOST_AUTO_TPL (RIGHTCAUCHYGREEN, transpose(deformationGradientTensor) * deformationGradientTensor);
BOOST_AUTO_TPL (firstInvariantC, trace( RIGHTCAUCHYGREEN ));
BOOST_AUTO_TPL (deformationGradientTensor, ( grad ( electrodETFESpace, *emDisp, 0) + value ( solid.activeMaterial()-> identity() ) ) );
BOOST_AUTO_TPL (RIGHTCAUCHYGREEN, transpose (deformationGradientTensor) * deformationGradientTensor);
BOOST_AUTO_TPL (firstInvariantC, trace ( RIGHTCAUCHYGREEN ) );
//BOOST_AUTO_TPL (fiber0, ( value( electrodETFESpace, *( monodomain -> fiberPtr() ) ) ));
#define fiber0 ( value( electrodETFESpace, *( monodomain -> fiberPtr() ) ) )
BOOST_AUTO_TPL (fiber, ( deformationGradientTensor * fiber0 ));
BOOST_AUTO_TPL (I4f, dot( fiber, fiber));
BOOST_AUTO_TPL (Ca, ( value( aETFESpace, *( monodomain -> globalSolution().at(3) ) ) ));
BOOST_AUTO_TPL (Ca2, ( Ca * Ca ));
BOOST_AUTO_TPL (dCa, ( Ca + value(-0.02155) ));
BOOST_AUTO_TPL (Gammaf, ( value( aETFESpace, *gammaf ) ));
BOOST_AUTO_TPL (GammaPlusOne, ( Gammaf + value(1.0) ));
BOOST_AUTO_TPL (Pa, ( value(-2.5) * eval(H, dCa ) * eval(H, dCa )/*eval( fl, I4f) */ * eval( flg, Gammaf) + value(1.0) / ( ( GammaPlusOne ) * ( GammaPlusOne ) * ( GammaPlusOne ) * ( GammaPlusOne ) ) ));
// #define Pa beta * eval( fl, I4f)
BOOST_AUTO_TPL (fiber, ( deformationGradientTensor * fiber0 ) );
BOOST_AUTO_TPL (I4f, dot ( fiber, fiber) );
BOOST_AUTO_TPL (Ca, ( value ( aETFESpace, * ( monodomain -> globalSolution().at (3) ) ) ) );
BOOST_AUTO_TPL (Ca2, ( Ca * Ca ) );
BOOST_AUTO_TPL (dCa, ( Ca + value (-0.02155) ) );
BOOST_AUTO_TPL (Gammaf, ( value ( aETFESpace, *gammaf ) ) );
BOOST_AUTO_TPL (GammaPlusOne, ( Gammaf + value (1.0) ) );
BOOST_AUTO_TPL (Pa, ( value (-2.5) * eval (H, dCa ) * eval (H, dCa ) /*eval( fl, I4f) */ * eval ( flg, Gammaf) + value (1.0) / ( ( GammaPlusOne ) * ( GammaPlusOne ) * ( GammaPlusOne ) * ( GammaPlusOne ) ) ) );
// #define Pa beta * eval( fl, I4f)
//#define dW1 ( ( I4f - value(1.0) ) )
BOOST_AUTO_TPL (dW1, ( ( Gammaf * Gammaf + 2.0 * Gammaf ) ));
BOOST_AUTO_TPL (dW, ( ( ( dW1 ) + value(1.0) / ( ( GammaPlusOne ) * ( GammaPlusOne ) * ( GammaPlusOne ) ) ) ));
BOOST_AUTO_TPL (dW1, ( ( Gammaf * Gammaf + 2.0 * Gammaf ) ) );
BOOST_AUTO_TPL (dW, ( ( ( dW1 ) + value (1.0) / ( ( GammaPlusOne ) * ( GammaPlusOne ) * ( GammaPlusOne ) ) ) ) );
// #define dgGammaf ( value(-1.0) + value(-2.0) / ( GammaPlusOne ) + value(2.0) * Gammaf * ( Gammaf + value(2.0) )* pow( GammaPlusOne, -3 ) )
// #define activationEquation value(-1.0) * ( Pa - ( value(2.0) * GammaPlusOne * firstInvariantC + dgGammaf * I4f ) * value( mu / 2.0 ) ) / beta
BOOST_AUTO_TPL (activationEquation, value(0.0005) * (Pa - dW) / ( Ca2 ));
BOOST_AUTO_TPL (activationEquation, value (0.0005) * (Pa - dW) / ( Ca2 ) );
vectorPtr_Type tmpRhsActivation ( new vector_Type ( rhsActivation -> map(), Repeated ) );
......
......@@ -530,16 +530,16 @@ int main (int argc, char** argv)
vectorPtr_Type solidFibers ( new vector_Type ( dFESpace -> map() ) );
ElectrophysiologyUtility::importFibers ( solidFibers,
solid_fiber_file,
localSolidMesh );
solid_fiber_file,
localSolidMesh );
solid.material() -> setFiberVector ( *solidFibers );
vectorPtr_Type solidSheet (new vector_Type ( dFESpace -> map() ) );
ElectrophysiologyUtility::importVectorField (solidSheet,
solid_sheet_file,
solid_sheet_field,
localSolidMesh );
solid_sheet_file,
solid_sheet_field,
localSolidMesh );
solid.material() -> setSheetVector (*solidSheet);
......
......@@ -212,7 +212,7 @@ int main ( int argc, char** argv )
std::string meshPath = parameterList.get ("mesh_path", "./");
meshPtr_Type meshPart (new mesh_Type ( Comm ) );
MeshUtility::loadMesh(meshPart, meshName, meshPath);
MeshUtility::loadMesh (meshPart, meshName, meshPath);
if (verbose)
{
......@@ -447,13 +447,13 @@ int main ( int argc, char** argv )
ExporterHDF5< RegionMesh <LinearTetra> > exporter;
exporter.setMeshProcId ( meshPart, Comm->MyPID() );
exporter.setPrefix ("rescalingGammaf");
exporter.setPostDir("./");
// exporter.addVariable ( ExporterData<mesh_Type>::ScalarField, "rescalingGammaf", uFESpace,
// solution, UInt (0) );
exporter.setPostDir ("./");
// exporter.addVariable ( ExporterData<mesh_Type>::ScalarField, "rescalingGammaf", uFESpace,
// solution, UInt (0) );
exporter.addVariable ( ExporterData<mesh_Type>::VectorField, "rescalingGammaf", uFESpace,
solution, UInt (0) );
exporter.postProcess ( 0 );
exporter.postProcess ( 0 );
if ( Comm->MyPID() == 0 )
{
......@@ -461,53 +461,53 @@ int main ( int argc, char** argv )
}
linearSolver.setRightHandSide (rhs);
linearSolver.solve (solution);
exporter.postProcess ( 1 );
// if ( Comm->MyPID() == 0 )
// {
// std::cout << "Transforming...\n";
// }
// Real p0 = 3.9981022451448112e-01;
// Real p1 = -1.7249368443147499e+00;
// Real p2 = 6.4695359113991486e+00;
// Real p3 = -1.9192450904013128e+01;
// Real p4 = 4.1297847160139170e+01;
// Real p5 = -5.9665612469298118e+01;
// Real p6 = 5.4040368339225651e+01;
// Real p7 = -2.7524476381745629e+01;
// Real p8 = 5.9999955063690678e+00;
//
// vector_Type tmp (solution -> map() );
// exporter.postProcess ( 0 );
// tmp *= 0.0;
// tmp += p0;
// tmp += ( p1 * *solution);
// tmp += ( p2 * *solution * *solution);
// tmp += ( p3 * *solution * *solution * *solution);
// tmp += ( p4 * *solution * *solution * *solution * *solution);
// tmp += ( p5 * *solution * *solution * *solution * *solution * *solution);
// tmp += ( p6 * *solution * *solution * *solution * *solution * *solution * *solution);
// tmp += ( p7 * *solution * *solution * *solution * *solution * *solution * *solution * *solution);
// tmp += ( p8 * *solution * *solution * *solution * *solution * *solution * *solution * *solution * *solution);
//
// *solution = tmp;
// exporter.postProcess ( 1 );
exporter.postProcess ( 1 );
// if ( Comm->MyPID() == 0 )
// {
// std::cout << "Transforming...\n";
// }
// Real p0 = 3.9981022451448112e-01;
// Real p1 = -1.7249368443147499e+00;
// Real p2 = 6.4695359113991486e+00;
// Real p3 = -1.9192450904013128e+01;
// Real p4 = 4.1297847160139170e+01;
// Real p5 = -5.9665612469298118e+01;
// Real p6 = 5.4040368339225651e+01;
// Real p7 = -2.7524476381745629e+01;
// Real p8 = 5.9999955063690678e+00;
//
// vector_Type tmp (solution -> map() );
// exporter.postProcess ( 0 );
// tmp *= 0.0;
// tmp += p0;
// tmp += ( p1 * *solution);
// tmp += ( p2 * *solution * *solution);
// tmp += ( p3 * *solution * *solution * *solution);
// tmp += ( p4 * *solution * *solution * *solution * *solution);
// tmp += ( p5 * *solution * *solution * *solution * *solution * *solution);
// tmp += ( p6 * *solution * *solution * *solution * *solution * *solution * *solution);
// tmp += ( p7 * *solution * *solution * *solution * *solution * *solution * *solution * *solution);
// tmp += ( p8 * *solution * *solution * *solution * *solution * *solution * *solution * *solution * *solution);
//
// *solution = tmp;
// exporter.postProcess ( 1 );
exporter.closeFile();
//
// vectorPtr_Type importedSolution ( new vector_Type ( solution -> map() ) );
//
// std::string filename = parameterList.get ("filename", "rescalingGammaf");
// std::string fieldname = parameterList.get ("fieldname", "rescalingGammaf");
// ElectrophysiologyUtility::importScalarField (importedSolution, filename, fieldname, meshPart);
//
// ExporterHDF5< RegionMesh <LinearTetra> > exporter2;
// exporter2.setMeshProcId ( meshPart, Comm->MyPID() );
// exporter2.setPrefix ("rescalingGammaf_refined1");
// exporter2.addVariable ( ExporterData<mesh_Type>::ScalarField, "rescalingGammaf", uFESpace,
// importedSolution, UInt (0) );
// exporter2.postProcess ( 0 );
// exporter2.closeFile();
//
// vectorPtr_Type importedSolution ( new vector_Type ( solution -> map() ) );
//
// std::string filename = parameterList.get ("filename", "rescalingGammaf");
// std::string fieldname = parameterList.get ("fieldname", "rescalingGammaf");
// ElectrophysiologyUtility::importScalarField (importedSolution, filename, fieldname, meshPart);
//
// ExporterHDF5< RegionMesh <LinearTetra> > exporter2;
// exporter2.setMeshProcId ( meshPart, Comm->MyPID() );
// exporter2.setPrefix ("rescalingGammaf_refined1");
// exporter2.addVariable ( ExporterData<mesh_Type>::ScalarField, "rescalingGammaf", uFESpace,
// importedSolution, UInt (0) );
// exporter2.postProcess ( 0 );
// exporter2.closeFile();
// ---------------------------------------------------------------
// We finalize the MPI session if MPI was used
// ---------------------------------------------------------------
......
......@@ -322,7 +322,7 @@ int main (int argc, char** argv)
std::cout << "\nsolid iterations: " << solidIter;
std::cout << "\nSolving coupled problem: " << coupling;
}
// Int k (0);
// Int k (0);
//// Real timeReac = 0.0;
//// Real timeDiff = 0.0;
......@@ -409,7 +409,7 @@ int main (int argc, char** argv)
expc->addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField, "disp", emSolverPtr -> activationPtr() -> FESpacePtrSolid(), emSolverPtr -> solidPtr() -> displacementPtr(), UInt (0) );
// UInt ntotalDof = emSolverPtr -> solidPtr() -> dispFESpace().dof().numTotalDof();
// UInt ntotalDof = emSolverPtr -> solidPtr() -> dispFESpace().dof().numTotalDof();
UInt nLocalDof = emSolverPtr -> solidPtr() -> displacement().epetraVector().MyLength();
UInt dim = nLocalDof / 3;
......@@ -467,104 +467,104 @@ int main (int argc, char** argv)
std::cout << "\n/////////////////////////////////////////////////////";
std::cout << "\n Time " << t;
std::cout << "\n/////////////////////////////////////////////////////";
}
if ( comm->MyPID() == 0 )
{
output << "\nnormres" << ", " << "normdisp" << ", time: " << t ;
}
Real normres = 2.0;
Real normdisp = 2.0;
iterk = 0;
while (normres > tol)
{
iterk++;
emSolverPtr -> solidPtr() -> displacement() = *Xoriginal;
emSolverPtr -> solidPtr() -> displacement() -= *Xk;
emSolverPtr -> solveSolid();
emSolverPtr -> exportSolution (t + 0.01 * iterk);
*res = *Xk;
*res += emSolverPtr -> solidPtr() -> displacement();
*res -= *Xoriginal;
normres = res -> normInf();
normdisp = emSolverPtr -> solidPtr() -> displacement().normInf();
//
if ( comm->MyPID() == 0 )
{
std::cout << "\n RESIDUAL : " << normres << ", Norm disp: " << normdisp ;
output << "\n" << normres << ", " << normdisp;
}
*Xk *= 0.0;
*Xk += *Xoriginal;
*Xk -= emSolverPtr -> solidPtr() -> displacement();
emSolverPtr -> solidPtr() -> mesh() -> meshTransformer().moveMesh ( *Xk, dim, 1);
*Xkp1 = emSolverPtr -> solidPtr() -> displacement();
//emSolverPtr -> solidPtr() -> mesh() -> meshTransformer().moveMesh( -emSolverPtr -> solidPtr() -> displacement(), dim, 0);
expc -> postProcess (t + 0.01 * iterk);
if (iterk > 15)
{
break;
}
}
}
output.close();
// //
// if ( comm->MyPID() == 0 )
// {
// cout << "\n---------------------------";
// cout << "\nInterpolating displacement ";
// cout << "\n---------------------------";
// }
// emSolverPtr -> updateMonodomain();
// }
//
//
// }
// emSolverPtr -> registerActivationTime(t, 0.8);
// matrixPtr_Type broydenMatrix( emSolverPtr -> solidPtr() -> );
// if( k % iter == 0 )
// {
// if ( comm->MyPID() == 0 )
// {
// cout << "\nExporting solutions";
// }
//emSolverPtr -> solidPtr() -> mesh() -> meshTransformer().moveMesh( -( emSolverPtr -> solidPtr() -> displacement()), 3);
//emSolverPtr -> exportSolution(2.0);
// }
//
// }
//
// if ( comm->MyPID() == 0 )
// {
// cout << "\nExporting Activation Time";
// }
iterk++;
emSolverPtr -> exportSolution (iterk);
emSolverPtr -> closeExporters();
expc -> closeFile();
std::cout << "\n\nThank you for using EMSolver.\nI hope to meet you again soon!\n All the best for your simulation :P\n " ;
// }
}
if ( comm->MyPID() == 0 )
{
output << "\nnormres" << ", " << "normdisp" << ", time: " << t ;
}
Real normres = 2.0;
Real normdisp = 2.0;
iterk = 0;
while (normres > tol)
{
iterk++;
emSolverPtr -> solidPtr() -> displacement() = *Xoriginal;
emSolverPtr -> solidPtr() -> displacement() -= *Xk;
emSolverPtr -> solveSolid();
emSolverPtr -> exportSolution (t + 0.01 * iterk);
*res = *Xk;
*res += emSolverPtr -> solidPtr() -> displacement();
*res -= *Xoriginal;
normres = res -> normInf();
normdisp = emSolverPtr -> solidPtr() -> displacement().normInf();
//
if ( comm->MyPID() == 0 )
{
std::cout << "\n RESIDUAL : " << normres << ", Norm disp: " << normdisp ;
output << "\n" << normres << ", " << normdisp;
}
*Xk *= 0.0;
*Xk += *Xoriginal;
*Xk -= emSolverPtr -> solidPtr() -> displacement();
emSolverPtr -> solidPtr() -> mesh() -> meshTransformer().moveMesh ( *Xk, dim, 1);
*Xkp1 = emSolverPtr -> solidPtr() -> displacement();
//emSolverPtr -> solidPtr() -> mesh() -> meshTransformer().moveMesh( -emSolverPtr -> solidPtr() -> displacement(), dim, 0);
expc -> postProcess (t + 0.01 * iterk);
if (iterk > 15)
{
break;
}
}
}
output.close();
// //
// if ( comm->MyPID() == 0 )
// {
// cout << "\n---------------------------";
// cout << "\nInterpolating displacement ";
// cout << "\n---------------------------";
// }
// emSolverPtr -> updateMonodomain();
// }
//
//
// }
// emSolverPtr -> registerActivationTime(t, 0.8);
// matrixPtr_Type broydenMatrix( emSolverPtr -> solidPtr() -> );
// if( k % iter == 0 )
// {
// if ( comm->MyPID() == 0 )
// {
// cout << "\nExporting solutions";
// }
//emSolverPtr -> solidPtr() -> mesh() -> meshTransformer().moveMesh( -( emSolverPtr -> solidPtr() -> displacement()), 3);
//emSolverPtr -> exportSolution(2.0);
// }
//
// }
//
// if ( comm->MyPID() == 0 )
// {
// cout << "\nExporting Activation Time";
// }
iterk++;
emSolverPtr -> exportSolution (iterk);
emSolverPtr -> closeExporters();
expc -> closeFile();
std::cout << "\n\nThank you for using EMSolver.\nI hope to meet you again soon!\n All the best for your simulation :P\n " ;
// }
#ifdef HAVE_MPI
MPI_Finalize();
MPI_Finalize();
#endif
return 0;
}
return 0;
}
......@@ -287,31 +287,31 @@ class orthonormalizeFibers
public:
typedef LifeV::VectorSmall<3> return_Type;
// return_Type operator() (const LifeV::VectorSmall<3>& v)
// {
// auto a(v);
// EMUtility::normalize(a);
// return a;
// }
// return_Type operator() (const LifeV::VectorSmall<3>& v)
// {
// auto a(v);
// EMUtility::normalize(a);
// return a;
// }
return_Type operator() (const LifeV::VectorSmall<3>& v)
{
auto a(v);
EMUtility::normalize(a, M_component);
return a;
auto a (v);
EMUtility::normalize (a, M_component);
return a;
}
return_Type operator() (const LifeV::VectorSmall<3>& v, const LifeV::VectorSmall<3>& w )
{
auto a(v);
auto b(w);
EMUtility::normalize(a);
EMUtility::orthonormalize(b, a);
return b;
auto a (v);
auto b (w);
EMUtility::normalize (a);
EMUtility::orthonormalize (b, a);
return b;
}
orthonormalizeFibers(LifeV::Int component = 0) : M_component(component) {}
~orthonormalizeFibers(){}
orthonormalizeFibers (LifeV::Int component = 0) : M_component (component) {}
~orthonormalizeFibers() {}
LifeV::Int M_component;
};
......@@ -322,13 +322,13 @@ class CrossProduct
public:
typedef LifeV::VectorSmall<3> return_Type;
return_Type operator()(const LifeV::VectorSmall<3>& v1, const LifeV::VectorSmall<3>& v2)
{
return Elasticity::crossProduct(v1, v2);
}
return_Type operator() (const LifeV::VectorSmall<3>& v1, const LifeV::VectorSmall<3>& v2)
{
return Elasticity::crossProduct (v1, v2);
}
CrossProduct() {}
~CrossProduct() {}
CrossProduct() {}
~CrossProduct() {}
};
......
This diff is collapsed.
This diff is collapsed.
......@@ -12,52 +12,54 @@
#include <lifev/core/array/MapEpetra.hpp>
namespace LifeV {
namespace LifeV
{
class ActiveStressActivation {
class ActiveStressActivation
{
public:
//! Distributed vector // For parallel usage
typedef VectorEpetra vector_Type;
typedef VectorEpetra vector_Type;
typedef boost::shared_ptr<VectorEpetra> vectorPtr_Type;
typedef boost::shared_ptr<VectorEpetra> vectorPtr_Type;
ActiveStressActivation(MapEpetra& map) : M_activationPtr( new vector_Type(map) ){}
ActiveStressActivation(const MapEpetra& map) : M_activationPtr( new vector_Type(map) ){}
ActiveStressActivation (MapEpetra& map) : M_activationPtr ( new vector_Type (map) ) {}
ActiveStressActivation (const MapEpetra& map) : M_activationPtr ( new vector_Type (map) ) {}
ActiveStressActivation(ActiveStressActivation& activation) : M_activationPtr( new vector_Type(activation.activation()) ){}
ActiveStressActivation (ActiveStressActivation& activation) : M_activationPtr ( new vector_Type (activation.activation() ) ) {}
inline ActiveStressActivation& operator=(ActiveStressActivation& activation)
{
M_activationPtr.reset( new vector_Type(activation.activation()) );
return *this;
}
inline ActiveStressActivation& operator= (ActiveStressActivation& activation)
{
M_activationPtr.reset ( new vector_Type (activation.activation() ) );
return *this;
}
virtual ~ActiveStressActivation() {}
virtual ~ActiveStressActivation() {}
inline VectorEpetra& activation()
{
return *M_activationPtr;
}
inline VectorEpetra& activation()