HeartSolver.hpp 2.41 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
#include <lifev/em/solver/EMSolver.hpp>
thomaskummer's avatar
thomaskummer committed
15
#include <lifev/em/solver/circulation/Circulation.hpp>
thomaskummer's avatar
thomaskummer committed
16
17


thomaskummer's avatar
thomaskummer committed
18
19
20
21

namespace LifeV
{

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

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

    
thomaskummer's avatar
thomaskummer committed
53
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
    
    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
97
98
99
100
};

}

thomaskummer's avatar
thomaskummer committed
101
#endif