fluid_problem_fem.h 6.8 KB
Newer Older
1
2
#ifndef FLUID_PROBLEM_FEM_H
#define FLUID_PROBLEM_FEM_H
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
3
#define D DIMENSION
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
4
#include <math.h>
5
#include "mesh.h"
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
6

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
7
#define  deps 1
8
9
10
11
12
13
14
15
16
struct FEElementStruct{
  int nlocal;
  int n[DIMENSION+1];
  void (*f)(const double *xi, double *f);
  void (*df)(const double *xi, const double dxidx[D][D], double f[][D]);
};

FEElement *fe_element_new(const char *etype);
void fe_element_free(FEElement *element);
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
17

Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
18
struct FEFieldsStruct{
19
20
21
22
  int local_size;
  int n;
  int nd[DIMENSION+1];
  FEElement **element;
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
23
};
24
25
26
27

FEFields *fe_fields_new(const char **etypes);
void fe_fields_f(const FEFields *fields, const double *xi, double *f);
void fe_fields_df(const FEFields *fields, const double *xi, const double dxdidx[D][D], double df[][D]);
Jonathan Lambrechts's avatar
wip    
Jonathan Lambrechts committed
28
29
30
31
void fe_fields_eval_grad_sf(const FEFields *fields, const Mesh *mesh, int iel, const double sf[], const double dsf[][DIMENSION], const double *v, double s[], double ds[][DIMENSION]);
void fe_fields_eval_sf(const FEFields *fields, const Mesh *mesh, int iel, const double sf[], const double *v, double s[]);
void fe_fields_eval_grad(const FEFields *fields, const double *xi, const double dxidx[D][D], const Mesh *mesh, int iel, const double *v, double s[], double ds[][DIMENSION]);
void fe_fields_eval(const FEFields *fields, const double *xi, const Mesh *mesh, int iel, const double *v, double s[]);
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
void fe_fields_add_to_local_vector(const FEFields *fields, const double *f0, const double f1[][D], const double *sf, const double dsf[][D], double jw, double *local_vector);
static void fe_fields_add_to_local_matrix(const FEFields *fields, const double *f00, const double f01[][D], const double f10[][D], const double f11[][D][D], const double *sf, const double dsf[][D], double jw, double *local_matrix) {
  int jc = 0;
  for (int jfield = 0; jfield < fields->n; ++jfield) {
    for (int jphi = 0; jphi < fields->element[jfield]->nlocal; ++jphi){
      int ic = 0;
      for (int ifield = 0; ifield < fields->n; ++ifield){
        for (int iphi = 0; iphi < fields->element[ifield]->nlocal; ++iphi){
          double m = 0;
          if (f00) {
            m += jw*sf[ic]*sf[jc]*f00[ifield*fields->n+jfield];
          }
          if (f10) {
            for (int id = 0; id < D; ++id) {
              m += jw*dsf[ic][id]*sf[jc]*f10[ifield*fields->n+jfield][id];
            }
          }
          if (f11) {
            for (int id = 0; id < D; ++id) {
              for (int jd = 0; jd < D; ++jd) {
                m += jw*dsf[ic][id]*dsf[jc][jd]*f11[ifield*fields->n+jfield][id][jd];
              }
            }
          }
          if (f01) {
            for (int jd = 0; jd < D; ++jd) {
              m += jw*sf[ic]*dsf[jc][jd]*f01[ifield*fields->n+jfield][jd];
            }
          }
          local_matrix[(ifield*3+iphi)*9+jfield*3+jphi] += m;
          ic ++;
        }
      }
      jc ++;
    }
  }
}
void fe_fields_free(FEFields *fields);

#if D==2
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
72
73
74
#define N_QP 3
#define N_LQP 2

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
75
static double LQP[][1] = {{0.21132486540518713}, {0.7886751345948129}};
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
76
77
78
static double LQW[] = {0.5,0.5};
static double QP[][2] = {{1./6,1./6},{1./6,2./3},{2./3,1./6}};
static double QW[] = {1./6,1./6,1./6};
79

Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
80
81
82
83
84
85
86
87
88
89
90
91
92
static double detDD(const double m[2][2])
{
  return m[0][0]*m[1][1] - m[0][1]*m[1][0];
}
static double invDD(const double m[2][2], double inv[2][2])
{
  double d = detDD(m);
  inv[0][0] = m[1][1]/d;
  inv[0][1] = -m[0][1]/d;
  inv[1][0] = -m[1][0]/d;
  inv[1][1] = m[0][0]/d;
  return d;
}
93

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
94
95
96
97
98
static const double invmassmatrix[3][3] = {
  {18,-6,-6},
  {-6,18,-6},
  {-6,-6,18}
};
Matthieu Constant's avatar
Matthieu Constant committed
99
100
101
102
103
static const double mass_matrix[3][3] = {
  {1./12,1./24,1./24},
  {1./24,1./12,1./24},
  {1./24,1./24,1./12}
};
104
105
106
107
108
109
110
111
112
113
114
115

static double get_normal(const Mesh *mesh, const int *interface,  double n[2]){
  const int *tri = &mesh->elements[interface[0]*3];
  const int *cl = elbnd2[interface[1]];
  double *x[2] = {&mesh->x[tri[cl[0]]*3], &mesh->x[tri[cl[1]]*3]};
  const double dx[2] = {x[1][0]-x[0][0], x[1][1]-x[0][1]};
  const double detbnd = sqrt(dx[0]*dx[0] + dx[1]*dx[1]);
  n[0] = dx[1]/detbnd;
  n[1] = -dx[0]/detbnd;
  return detbnd;
}

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
116
117
118
119
120
121
122
static void get_normal_and_det(double *x[2], double n[2], double *det){
  const double dx[2] = {x[1][0]-x[0][0], x[1][1]-x[0][1]};
  const double detbnd = sqrt(dx[0]*dx[0] + dx[1]*dx[1]);
  n[0] = dx[1]/detbnd;
  n[1] = -dx[0]/detbnd;
  *det = detbnd;
}
123

Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
#else

#define N_QP 5
#define N_LQP 3
static double QP[][3] = {
  {0.25, 0.25, 0.25},
  {0.5, 1./6, 1./6},
  {1./6, 0.5, 1./6},
  {1./6, 1./6, 0.5},
  {1./6, 1./6, 1./6}
};
static double LQP[][2] = {{1./6,1./6},{1./6,2./3},{2./3,1./6}};
static double QW[] = {-16./120, 9./120, 9./120, 9./120, 9./120};
static double LQW[] = {1./6,1./6,1./6};
static double detDD(const double m[3][3])
{
  return m[0][0]*(m[1][1]*m[2][2]-m[2][1]*m[1][2])
        -m[1][0]*(m[0][1]*m[2][2]-m[2][1]*m[0][2])
        +m[2][0]*(m[0][1]*m[1][2]-m[1][1]*m[0][2]);
}
static double invDD(const double m[3][3], double inv[3][3]) {
  double d = detDD(m);
  inv[0][0] = (m[1][1]*m[2][2]-m[2][1]*m[1][2])/d;
  inv[0][1] = -(m[0][1]*m[2][2]-m[2][1]*m[0][2])/d;
  inv[0][2] = (m[0][1]*m[1][2]-m[1][1]*m[0][2])/d;
  inv[1][0] = -(m[1][0]*m[2][2]-m[2][0]*m[1][2])/d;
  inv[1][1] = (m[0][0]*m[2][2]-m[2][0]*m[0][2])/d;
  inv[1][2] = -(m[0][0]*m[1][2]-m[1][0]*m[0][2])/d;
  inv[2][0] = (m[1][0]*m[2][1]-m[2][0]*m[1][1])/d;
  inv[2][1] = -(m[0][0]*m[2][1]-m[2][0]*m[0][1])/d;
  inv[2][2] = (m[0][0]*m[1][1]-m[1][0]*m[0][1])/d;
  return d;
};

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
158
159
160
161
162
163
static const double invmassmatrix[4][4] = {
  { 96,-24,-24,-24},
  {-24, 96,-24,-24},
  {-24,-24, 96,-24},
  {-24,-24,-24, 96}
};
Matthieu Constant's avatar
Matthieu Constant committed
164
165
166
167
168
169
static const double mass_matrix[4][4] = {
  {1./60 ,1./120,1./120,1./120},
  {1./120,1./60 ,1./120,1./120},
  {1./120,1./120,1./60 ,1./120},
  {1./120,1./120,1./120,1./60 }
};
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
170
171
172
173
174
175
176
177
178
179
180
181
static void get_normal_and_det(double *x[3], double n[3], double *det) {
  const double dx[2][3] = {{x[1][0]-x[0][0], x[1][1]-x[0][1], x[1][2]-x[0][2]},{x[2][0]-x[0][0], x[2][1]-x[0][1], x[2][2]-x[0][2]}};
  const double nn[3] = {
    dx[0][1]*dx[1][2] - dx[0][2]*dx[1][1],
    dx[0][2]*dx[1][0] - dx[0][0]*dx[1][2],
    dx[0][0]*dx[1][1] - dx[0][1]*dx[1][0]};
  const double detbnd = sqrt(nn[0]*nn[0] + nn[1]*nn[1] + nn[2]*nn[2]);
  n[0] = -nn[0]/detbnd;
  n[1] = -nn[1]/detbnd;
  n[2] = -nn[2]/detbnd;
  *det = detbnd;
}
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
182

183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
static double get_normal(const Mesh *mesh, const int *interface,  double n[3]){
  const int *tet = &mesh->elements[interface[0]*4];
  const int *cl = elbnd2[interface[1]];
  double *x[3] = {&mesh->x[tet[cl[0]]*3], &mesh->x[tet[cl[1]]*3],  &mesh->x[tet[cl[2]]*3]};
  const double dx[2][3] = {{x[1][0]-x[0][0], x[1][1]-x[0][1], x[1][2]-x[0][2]},{x[2][0]-x[0][0], x[2][1]-x[0][1], x[2][2]-x[0][2]}};
  const double nn[3] = {
    dx[0][1]*dx[1][2] - dx[0][2]*dx[1][1],
    dx[0][2]*dx[1][0] - dx[0][0]*dx[1][2],
    dx[0][0]*dx[1][1] - dx[0][1]*dx[1][0]};
  const double detbnd = sqrt(nn[0]*nn[0] + nn[1]*nn[1] + nn[2]*nn[2]);
  n[0] = -nn[0]/detbnd;
  n[1] = -nn[1]/detbnd;
  n[2] = -nn[2]/detbnd;
  return detbnd;
}
Jonathan Lambrechts's avatar
broken    
Jonathan Lambrechts committed
198

199
200
#endif
#endif