HeartSolver.hpp 2.38 KB
Newer Older
thomaskummer's avatar
thomaskummer committed
1
2
3
4
5
6
7
8
9
10
11
12
13
//
//  HeartSolver.cpp
//  LifeV
//
//  Created by Thomas Kummer on 03.05.16.
//  Copyright © 2016 Thomas Kummer. All rights reserved.
//

#ifndef _HEARTDATA_H_
#define _HEARTDATA_H_


#include <stdio.h>
thomaskummer's avatar
thomaskummer committed
14
15
16
#include <lifev/em/solver/EMSolver.hpp>


thomaskummer's avatar
thomaskummer committed
17
18
19
20

namespace LifeV
{

thomaskummer's avatar
thomaskummer committed
21
template<class meshType>
thomaskummer's avatar
thomaskummer committed
22
23
24
25
class HeartSolver {
   
public:
    
thomaskummer's avatar
thomaskummer committed
26
27
    
    HeartSolver(EMSolver<meshType,EMMonodomainSolver<meshType> > & emSolver,  CirculationSolver & circulationSolver):
thomaskummer's avatar
thomaskummer committed
28
        M_emSolver          (emSolver),
thomaskummer's avatar
thomaskummer committed
29
30
31
32
33
        M_circulationSolver (circulationSolver)
    {};
    
    virtual ~HeartSolver() {};
    
thomaskummer's avatar
thomaskummer committed
34
35

    
thomaskummer's avatar
thomaskummer committed
36
37
38
39
40
    
    
    
    
    
thomaskummer's avatar
thomaskummer committed
41
protected:
thomaskummer's avatar
thomaskummer committed
42
    
thomaskummer's avatar
thomaskummer committed
43
    
thomaskummer's avatar
thomaskummer committed
44
    EMSolver<meshType,EMMonodomainSolver<meshType> > M_emSolver;
thomaskummer's avatar
thomaskummer committed
45
46
47
48
49
    CirculationSolver M_circulationSolver;
    
    
    
    
thomaskummer's avatar
thomaskummer committed
50
51
52
53
    VectorSmall<2> M_pressure;
    VectorSmall<2> M_volume;

    
thomaskummer's avatar
thomaskummer committed
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
    
    Real
    externalPower ( const VectorEpetra& dispCurrent,
                    const VectorEpetra& dispPrevious,
                    const boost::shared_ptr<ETFESpace< RegionMesh<LinearTetra>, MapEpetra, 3, 3 > > dispETFESpace,
                    Real pressure,
                    Real dt,
                    const unsigned int bdFlag) const
    {
        VectorEpetra traction ( dispCurrent.map() );
        VectorEpetra velocity ( (dispCurrent - dispPrevious) / dt );
        
        MatrixSmall<3,3> Id;
        Id(0,0) = 1.; Id(0,1) = 0., Id(0,2) = 0.;
        Id(1,0) = 0.; Id(1,1) = 1., Id(1,2) = 0.;
        Id(2,0) = 0.; Id(2,1) = 0., Id(2,2) = 1.;
        
        {
            using namespace ExpressionAssembly;
            
            auto I = value(Id);
            auto Grad_u = grad( dispETFESpace, dispCurrent, 0);
            auto F =  Grad_u + I;
            auto FmT = minusT(F);
            auto J = det(F);
            auto p = value(pressure);
            
            QuadratureBoundary myBDQR (buildTetraBDQR (quadRuleTria7pt) );
            
            integrate ( boundary ( dispETFESpace->mesh(), bdFlag),
                       myBDQR,
                       dispETFESpace,
                       p * J * dot( FmT * Nface,  phi_i)
                       //p * J * dot( FmT * Nface,  phi_i)
                       //value(-1.0) * J * dot (vE1, FmT * Nface) * phi_i) >> intergral
                       ) >> traction;
            
            traction.globalAssemble();
        }
        
        return traction.dot(velocity);
    }
    
    
thomaskummer's avatar
thomaskummer committed
98
99
100
101
};

}

thomaskummer's avatar
thomaskummer committed
102
#endif