Trilinos based (stochastic) FEM solvers
fepp.cpp
Go to the documentation of this file.
1 /*
2 Brian Staber (brian.staber@gmail.com)
3 */
4 
5 #include "fepp.hpp"
6 
7 void tri3::shape_functions(Epetra_SerialDenseVector & N, double & xi, double & eta){
8  N(0) = 1.0 - xi - eta;
9  N(1) = xi;
10  N(2) = eta;
11 }
12 void tri3::d_shape_functions(Epetra_SerialDenseMatrix & D, double & xi, double & eta){
13  D(0,0) = -1.0; D(0,1) = -1.0;
14  D(1,0) = 1.0; D(1,1) = 0.0;
15  D(2,0) = 0.0; D(2,1) = 1.0;
16 }
17 void tri3::dX_shape_functions(Epetra_SerialDenseMatrix & D, Epetra_SerialDenseMatrix & DX, double & jac, Epetra_SerialDenseMatrix & JacobianMatrix){
18 
19  Epetra_SerialDenseMatrix InverseJacobianMatrix(2,2);
20  InverseJacobianMatrix(0,0) = (1.0/jac)*JacobianMatrix(1,1);
21  InverseJacobianMatrix(1,1) = (1.0/jac)*JacobianMatrix(0,0);
22  InverseJacobianMatrix(0,1) = -(1.0/jac)*JacobianMatrix(0,1);
23  InverseJacobianMatrix(1,0) = -(1.0/jac)*JacobianMatrix(1,0);
24 
25  DX.Multiply('N','N',1.0,D,InverseJacobianMatrix,0.0);
26 }
27 
28 void quad4::shape_functions(Epetra_SerialDenseVector & N, double & xi, double & eta){
29  N(0) = (1.0/4.0)*(1.0-xi)*(1.0-eta);
30  N(1) = (1.0/4.0)*(1.0+xi)*(1.0-eta);
31  N(2) = (1.0/4.0)*(1.0+xi)*(1.0+eta);
32  N(3) = (1.0/4.0)*(1.0-xi)*(1.0+eta);
33 }
34 void quad4::d_shape_functions(Epetra_SerialDenseMatrix & D, double & xi, double & eta){
35  D(0,0) = -(1.0/4.0)*(1.0-eta); D(0,1) = -(1.0/4.0)*(1.0-xi);
36  D(1,0) = (1.0/4.0)*(1.0-eta); D(1,1) = -(1.0/4.0)*(1.0+xi);
37  D(2,0) = (1.0/4.0)*(1.0+eta); D(2,1) = (1.0/4.0)*(1.0+xi);
38  D(3,0) = -(1.0/4.0)*(1.0+eta); D(3,1) = (1.0/4.0)*(1.0-xi);
39 }
40 void quad4::dX_shape_functions(Epetra_SerialDenseMatrix & D, Epetra_SerialDenseMatrix & DX, double & jac, Epetra_SerialDenseMatrix & JacobianMatrix){
41 
42  Epetra_SerialDenseMatrix InverseJacobianMatrix(2,2);
43  InverseJacobianMatrix(0,0) = (1.0/jac)*JacobianMatrix(1,1);
44  InverseJacobianMatrix(1,1) = (1.0/jac)*JacobianMatrix(0,0);
45  InverseJacobianMatrix(0,1) = -(1.0/jac)*JacobianMatrix(0,1);
46  InverseJacobianMatrix(1,0) = -(1.0/jac)*JacobianMatrix(1,0);
47 
48  DX.Multiply('N','N',1.0,D,InverseJacobianMatrix,0.0);
49 }
50 
51 void tri6::shape_functions(Epetra_SerialDenseVector & N, double & xi, double & eta){
52  double lambda = 1.0 - xi - eta;
53  N(0) = lambda*(2.0*lambda-1.0);
54  N(1) = xi*(2.0*xi-1.0);
55  N(2) = eta*(2.0*eta-1.0);
56  N(3) = 4.0*lambda*xi;
57  N(4) = 4.0*eta*xi;
58  N(5) = 4.0*lambda*eta;
59 }
60 
61 void tri6::d_shape_functions(Epetra_SerialDenseMatrix & D, double & xi, double & eta){
62  D(0,0) = 4.0*xi + 4.0*eta - 3.0;
63  D(1,0) = 4.0*xi - 1.0;
64  D(2,0) = 0.0;
65  D(3,0) = 4.0 - 8.0*xi - 4.0*eta;
66  D(4,0) = 4.0*eta;
67  D(5,0) = -4.0*eta;
68 
69  D(0,1) = 4.0*eta + 4.0*xi - 3.0;
70  D(1,1) = 0.0;
71  D(2,1) = 4.0*eta - 1.0;
72  D(3,1) = -4.0*xi;
73  D(4,1) = 4.0*xi;
74  D(5,1) = 4.0 - 4.0*xi - 8.0*eta;
75 }
76 
77 void tri6::dX_shape_functions(Epetra_SerialDenseMatrix & D, Epetra_SerialDenseMatrix & DX, double & jac, Epetra_SerialDenseMatrix & JacobianMatrix){
78 
79  Epetra_SerialDenseMatrix InverseJacobianMatrix(2,2);
80  InverseJacobianMatrix(0,0) = (1.0/jac)*JacobianMatrix(1,1);
81  InverseJacobianMatrix(1,1) = (1.0/jac)*JacobianMatrix(0,0);
82  InverseJacobianMatrix(0,1) = -(1.0/jac)*JacobianMatrix(0,1);
83  InverseJacobianMatrix(1,0) = -(1.0/jac)*JacobianMatrix(1,0);
84 
85  DX.Multiply('N','N',1.0,D,InverseJacobianMatrix,0.0);
86 }
87 
88 void hexa8::shape_functions(Epetra_SerialDenseVector & N, double & xi, double & eta, double & zeta){
89  N(0) = (1.0/8.0)*(1.0-xi)*(1.0-eta)*(1.0-zeta);
90  N(1) = (1.0/8.0)*(1.0+xi)*(1.0-eta)*(1.0-zeta);
91  N(2) = (1.0/8.0)*(1.0+xi)*(1.0+eta)*(1.0-zeta);
92  N(3) = (1.0/8.0)*(1.0-xi)*(1.0+eta)*(1.0-zeta);
93  N(4) = (1.0/8.0)*(1.0-xi)*(1.0-eta)*(1.0+zeta);
94  N(5) = (1.0/8.0)*(1.0+xi)*(1.0-eta)*(1.0+zeta);
95  N(6) = (1.0/8.0)*(1.0+xi)*(1.0+eta)*(1.0+zeta);
96  N(7) = (1.0/8.0)*(1.0-xi)*(1.0+eta)*(1.0+zeta);
97 }
98 
99 void hexa8::d_shape_functions(Epetra_SerialDenseMatrix & D, double & xi, double & eta, double & zeta){
100  D(0,0) = -(1.0/8.0)*(1.0-eta)*(1.0-zeta); D(0,1) = -(1.0/8.0)*(1.0-xi)*(1.0-zeta); D(0,2) = -(1.0/8.0)*(1.0-xi)*(1.0-eta);
101  D(1,0) = (1.0/8.0)*(1.0-eta)*(1.0-zeta); D(1,1) = -(1.0/8.0)*(1.0+xi)*(1.0-zeta); D(1,2) = -(1.0/8.0)*(1.0+xi)*(1.0-eta);
102  D(2,0) = (1.0/8.0)*(1.0+eta)*(1.0-zeta); D(2,1) = (1.0/8.0)*(1.0+xi)*(1.0-zeta); D(2,2) = -(1.0/8.0)*(1.0+xi)*(1.0+eta);
103  D(3,0) = -(1.0/8.0)*(1.0+eta)*(1.0-zeta); D(3,1) = (1.0/8.0)*(1.0-xi)*(1.0-zeta); D(3,2) = -(1.0/8.0)*(1.0-xi)*(1.0+eta);
104  D(4,0) = -(1.0/8.0)*(1.0-eta)*(1.0+zeta); D(4,1) = -(1.0/8.0)*(1.0-xi)*(1.0+zeta); D(4,2) = (1.0/8.0)*(1.0-xi)*(1.0-eta);
105  D(5,0) = (1.0/8.0)*(1.0-eta)*(1.0+zeta); D(5,1) = -(1.0/8.0)*(1.0+xi)*(1.0+zeta); D(5,2) = (1.0/8.0)*(1.0+xi)*(1.0-eta);
106  D(6,0) = (1.0/8.0)*(1.0+eta)*(1.0+zeta); D(6,1) = (1.0/8.0)*(1.0+xi)*(1.0+zeta); D(6,2) = (1.0/8.0)*(1.0+xi)*(1.0+eta);
107  D(7,0) = -(1.0/8.0)*(1.0+eta)*(1.0+zeta); D(7,1) = (1.0/8.0)*(1.0-xi)*(1.0+zeta); D(7,2) = (1.0/8.0)*(1.0-xi)*(1.0+eta);
108 }
109 
110 void tetra4::shape_functions(Epetra_SerialDenseVector & N, double & xi, double & eta, double & zeta){
111  double lambda = 1 - xi - eta - zeta;
112  N(0) = lambda;
113  N(1) = xi;
114  N(2) = eta;
115  N(3) = zeta;
116 }
117 
118 void tetra4::d_shape_functions(Epetra_SerialDenseMatrix & D, double & xi, double & eta, double & zeta){
119  D(0,0) = -1.0; D(0,1) = -1.0; D(0,2) = -1.0;
120  D(1,0) = 1.0; D(1,1) = 0.0; D(1,2) = 0.0;
121  D(2,0) = 0.0; D(2,1) = 1.0; D(2,2) = 0.0;
122  D(3,0) = 0.0; D(3,1) = 0.0; D(3,2) = 1.0;
123 }
124 
125 void tetra10::shape_functions(Epetra_SerialDenseVector & N, double & xi, double & eta, double & zeta){
126  double lambda = 1 - xi - eta - zeta;
127  N(0) = lambda*(2.0*lambda-1.0);
128  N(1) = xi*(2.0*xi-1.0);
129  N(2) = eta*(2.0*eta-1.0);
130  N(3) = zeta*(2.0*zeta-1.0);
131  N(4) = 4*xi*lambda;
132  N(5) = 4*xi*eta;
133  N(6) = 4*eta*lambda;
134  N(7) = 4*zeta*lambda;
135  N(8) = 4*eta*zeta;
136  N(9) = 4*xi*zeta;
137 }
138 
139 void tetra10::d_shape_functions(Epetra_SerialDenseMatrix & D, double & xi, double & eta, double & zeta){
140  D(0,0) = 4.0*eta + 4.0*xi + 4.0*zeta - 3.0; D(0,1) = 4.0*eta + 4.0*xi + 4.0*zeta - 3.0; D(0,2) = 4.0*eta + 4.0*xi + 4.0*zeta - 3.0;
141  D(1,0) = 4.0*xi - 1.0; D(1,1) = 0.0; D(1,2) = 0.0;
142  D(2,0) = 0.0; D(2,1) = 4.0*eta - 1.0; D(2,2) = 0.0;
143  D(3,0) = 0.0; D(3,1) = 0.0; D(3,2) = 4.0*zeta - 1.0;
144  D(4,0) = 4.0 - 8.0*xi - 4.0*zeta - 4.0*eta; D(4,1) = -4.0*xi; D(4,2) = -4.0*xi;
145  D(5,0) = 4.0*eta; D(5,1) = 4.0*xi; D(5,2) = 0.0;
146  D(6,0) = -4.0*eta; D(6,1) = 4.0 - 4.0*xi - 4.0*zeta - 8.0*eta; D(6,2) = -4.0*eta;
147  D(7,0) = -4.0*zeta; D(7,1) = -4.0*zeta; D(7,2) = 4.0 - 4.0*xi - 8.0*zeta - 4.0*eta;
148  D(8,0) = 0.0; D(8,1) = 4.0*zeta; D(8,2) = 4.0*eta;
149  D(9,0) = 4.0*zeta; D(9,1) = 0.0; D(9,2) = 4.0*xi;
150 }
151 
152 void dX_shape_functions(Epetra_SerialDenseMatrix & D, Epetra_SerialDenseMatrix JacobianMatrix, double & jac, Epetra_SerialDenseMatrix & DX)
153 {
154  Epetra_SerialDenseMatrix InverseJacobianMatrix(3,3);
155 
156  InverseJacobianMatrix(0,0) = (1.0/jac)*(JacobianMatrix(1,1)*JacobianMatrix(2,2)-JacobianMatrix(1,2)*JacobianMatrix(2,1));
157  InverseJacobianMatrix(0,1) = (1.0/jac)*(JacobianMatrix(0,2)*JacobianMatrix(2,1)-JacobianMatrix(0,1)*JacobianMatrix(2,2));
158  InverseJacobianMatrix(0,2) = (1.0/jac)*(JacobianMatrix(0,1)*JacobianMatrix(1,2)-JacobianMatrix(0,2)*JacobianMatrix(1,1));
159 
160  InverseJacobianMatrix(1,0) = (1.0/jac)*(JacobianMatrix(1,2)*JacobianMatrix(2,0)-JacobianMatrix(1,0)*JacobianMatrix(2,2));
161  InverseJacobianMatrix(1,1) = (1.0/jac)*(JacobianMatrix(0,0)*JacobianMatrix(2,2)-JacobianMatrix(0,2)*JacobianMatrix(2,0));
162  InverseJacobianMatrix(1,2) = (1.0/jac)*(JacobianMatrix(0,2)*JacobianMatrix(1,0)-JacobianMatrix(0,0)*JacobianMatrix(1,2));
163 
164  InverseJacobianMatrix(2,0) = (1.0/jac)*(JacobianMatrix(1,0)*JacobianMatrix(2,1)-JacobianMatrix(1,1)*JacobianMatrix(2,0));
165  InverseJacobianMatrix(2,1) = (1.0/jac)*(JacobianMatrix(0,1)*JacobianMatrix(2,0)-JacobianMatrix(0,0)*JacobianMatrix(2,1));
166  InverseJacobianMatrix(2,2) = (1.0/jac)*(JacobianMatrix(0,0)*JacobianMatrix(1,1)-JacobianMatrix(0,1)*JacobianMatrix(1,0));
167 
168  DX.Multiply('N','N',1.0,D,InverseJacobianMatrix,0.0);
169 }
170 
171 void jacobian_matrix(Epetra_SerialDenseMatrix & X, Epetra_SerialDenseMatrix & D, Epetra_SerialDenseMatrix & JacobianMatrix){
172 JacobianMatrix.Multiply('N','N',1.0,X,D,0.0);
173 }
174 
175 void jacobian_det(Epetra_SerialDenseMatrix & JacobianMatrix, double & jac){
176  jac = fabs(JacobianMatrix(0,0)*JacobianMatrix(1,1)*JacobianMatrix(2,2)-JacobianMatrix(0,0)*JacobianMatrix(1,2)*JacobianMatrix(2,1)-JacobianMatrix(0,1)*JacobianMatrix(1,0)*JacobianMatrix(2,2)+JacobianMatrix(0,1)*JacobianMatrix(1,2)*JacobianMatrix(2,0)+JacobianMatrix(0,2)*JacobianMatrix(1,0)*JacobianMatrix(2,1)-JacobianMatrix(0,2)*JacobianMatrix(1,1)*JacobianMatrix(2,0) );
177 }
178 
179 void jacobian_det_faces(Epetra_SerialDenseMatrix & JacobianMatrix, double & jac){
180  jac = fabs(JacobianMatrix(0,0)*JacobianMatrix(1,1) - JacobianMatrix(0,1)*JacobianMatrix(0,1));
181 }
182 
183 void gauss_points_tri1(Epetra_SerialDenseVector & weight, Epetra_SerialDenseVector & xi, Epetra_SerialDenseVector & eta){
184  weight.Resize(1); xi.Resize(1); eta.Resize(1);
185  weight(0) = 1.0/2.0;
186  xi(0) = 1.0/3.0;
187  eta(0) = 1.0/3.0;
188 }
189 void gauss_points_tri3(Epetra_SerialDenseVector & weight, Epetra_SerialDenseVector & xi, Epetra_SerialDenseVector & eta){
190  weight.Resize(3); xi.Resize(3); eta.Resize(3);
191  weight(0) = 1.0/6.0; weight(1) = 1.0/6.0; weight(2) = 1.0/6.0;
192  xi(0) = 1.0/6.0; xi(1) = 2.0/3.0; xi(2) = 1.0/6.0;
193  eta(0) = 1.0/6.0; eta(1) = 1.0/6.0; eta(2) = 2.0/3.0;
194 }
195 void gauss_points_tri4(Epetra_SerialDenseVector & weight, Epetra_SerialDenseVector & xi, Epetra_SerialDenseVector & eta){
196  weight.Resize(4); xi.Resize(4); eta.Resize(4);
197  weight(0) = -27.0/96.0; weight(1) = 25.0/96.0; weight(2) = 25.0/96.0; weight(3) = 25.0/96.0;
198  xi(0) = 1.0/3.0; xi(1) = 1.0/5.0; xi(2) = 3.0/5.0; xi(3) = 1.0/5.0;
199  eta(0) = 1.0/3.0; eta(1) = 1.0/5.0; eta(2) = 1.0/5.0; eta(3) = 3.0/5.0;
200 }
201 
202 void gauss_points_quad1(Epetra_SerialDenseVector & weight, Epetra_SerialDenseVector & xi, Epetra_SerialDenseVector & eta){
203  weight.Resize(1); xi.Resize(1); eta.Resize(1);
204  weight(0) = 4.0;
205  xi(0) = 0.0;
206  eta(0) = 0.0;
207 }
208 void gauss_points_quad4(Epetra_SerialDenseVector & weight, Epetra_SerialDenseVector & xi, Epetra_SerialDenseVector & eta){
209  weight.Resize(4); xi.Resize(4); eta.Resize(4);
210  double a = 1.0/std::sqrt(3.0);
211  weight(0) = 1.0; weight(1) = 1.0; weight(2) = 1.0; weight(3) = 1.0;
212  xi(0) = -a; xi(1) = a; xi(2) = -a; xi(3) = a;
213  eta(0) = -a; eta(1) = -a; eta(2) = a; eta(3) = a;
214 }
215 void gauss_points_quad9(Epetra_SerialDenseVector & weight, Epetra_SerialDenseVector & xi, Epetra_SerialDenseVector & eta){
216  weight.Resize(9); xi.Resize(9); eta.Resize(9);
217  double a = std::sqrt(3.0/5.0);
218  weight(0) = 25.0/81.0; weight(1) = 40.0/81.0; weight(2) = 25.0/81.0; weight(3) = 40.0/81.0; weight(4) = 64.0/81.0; weight(5) = 40.0/81.0; weight(6) = 25.0/81.0; weight(7) = 40.0/81.0; weight(8) = 25.0/81.0;
219  xi(0) = -a; xi(1) = 0.0; xi(2) = a; xi(3) = -a; xi(4) = 0.0; xi(5) = a; xi(6) = -a; xi(7) = 0.0; xi(8) = a;
220  eta(0) = -a; eta(1) = -a; eta(2) = -a; eta(3) = 0.0; eta(4) = 0.0; eta(5) = 0.0; eta(6) = a; eta(7) = a; eta(8) = a;
221 }
222 
223 void gauss_points_hexa4(Epetra_SerialDenseVector & weight, Epetra_SerialDenseVector & xi, Epetra_SerialDenseVector & eta, Epetra_SerialDenseVector & zeta){
224  weight.Resize(4); xi.Resize(4); eta.Resize(4); zeta.Resize(4);
225  double alpha = std::sqrt(2.0/3.0);
226  double beta = std::sqrt(1.0/3.0);
227  weight(0) = 2.0; weight(1) = 2.0; weight(2) = 2.0; weight(3) = 2.0;
228  xi(0) = 0.0; eta(0) = alpha; zeta(0) = -beta;
229  xi(1) = 0.0; eta(1) = -alpha; zeta(1) = -beta;
230  xi(2) = alpha; eta(2) = 0.0; zeta(2) = beta;
231  xi(3)= -alpha; eta(3) = 0.0; zeta(3) = beta;
232 }
233 
234 void gauss_points_hexa8(Epetra_SerialDenseVector & weight, Epetra_SerialDenseVector & xi, Epetra_SerialDenseVector & eta, Epetra_SerialDenseVector & zeta){
235  weight.Resize(8); xi.Resize(8); eta.Resize(8); zeta.Resize(8);
236 
237  double a = 1.0/std::sqrt(3.0);
238  xi(0)=-a; eta(0)=-a; zeta(0)=-a; weight(0)=1.0;
239  xi(1)=-a; eta(1)=-a; zeta(1)=a; weight(1)=1.0;
240  xi(2)=-a; eta(2)=a; zeta(2)=-a; weight(2)=1.0;
241  xi(3)=-a; eta(3)=a; zeta(3)=a; weight(3)=1.0;
242  xi(4)= a; eta(4)=-a; zeta(4)=-a; weight(4)=1.0;
243  xi(5)= a; eta(5)=-a; zeta(5)=a; weight(5)=1.0;
244  xi(6)= a; eta(6)=a; zeta(6)=-a; weight(6)=1.0;
245  xi(7)= a; eta(7)=a; zeta(7)=a; weight(7)=1.0;
246 }
247 
248 void gauss_points_hexa27(Epetra_SerialDenseVector & weight, Epetra_SerialDenseVector & xi, Epetra_SerialDenseVector & eta, Epetra_SerialDenseVector & zeta){
249  weight.Resize(27); xi.Resize(27); eta.Resize(27); zeta.Resize(27);
250 
251  double a = std::sqrt(3.0/5.0); double c1 = 5.0/9.0; double c2 = 8.0/9.0;
252  xi(0) = -a; eta(0) = -a; zeta(0) = -a; weight(0) = c1*c1*c1;
253  xi(1) = -a; eta(1) = -a; zeta(1) = 0.0; weight(1) = c1*c1*c2;
254  xi(2) = -a; eta(2) = -a; zeta(2) = a; weight(2) = c1*c1*c1;
255  xi(3) = -a; eta(3) = 0.0; zeta(3) = -a; weight(3) = c1*c1*c2;
256  xi(4) = -a; eta(4) = 0.0; zeta(4) = 0.0; weight(4) = c1*c2*c2;
257  xi(5) = -a; eta(5) = 0.0; zeta(5) = a; weight(5) = c1*c1*c2;
258  xi(6) = -a; eta(6) = a; zeta(6) = -a; weight(6) = c1*c1*c1;
259  xi(7) = -a; eta(7) = a; zeta(7) = 0.0; weight(7) = c1*c1*c2;
260  xi(8) = -a; eta(8) = a; zeta(8) = a; weight(8) = c1*c1*c1;
261  xi(9) = 0.0; eta(9) = -a; zeta(9) = -a; weight(9) = c1*c1*c2;
262  xi(10) = 0.0; eta(10) = -a; zeta(10) = 0.0; weight(10) = c1*c2*c2;
263  xi(11) = 0.0; eta(11) = -a; zeta(11) = a; weight(11) = c1*c1*c2;
264  xi(12) = 0.0; eta(12) = 0.0; zeta(12) = -a; weight(12) = c1*c2*c2;
265  xi(13) = 0.0; eta(13) = 0.0; zeta(13) = 0.0; weight(13) = c2*c2*c2;
266  xi(14) = 0.0; eta(14) = 0.0; zeta(14) = a; weight(14) = c1*c2*c2;
267  xi(15) = 0.0; eta(15) = a; zeta(15) = -a; weight(15) = c1*c1*c2;
268  xi(16) = 0.0; eta(16) = a; zeta(16) = 0.0; weight(16) = c1*c2*c2;
269  xi(17) = 0.0; eta(17) = a; zeta(17) = a; weight(17) = c1*c1*c2;
270  xi(18) = a; eta(18) = -a; zeta(18) = -a; weight(18) = c1*c1*c1;
271  xi(19) = a; eta(19) = -a; zeta(19) = 0.0; weight(19) = c1*c1*c2;
272  xi(20) = a; eta(20) = -a; zeta(20) = a; weight(20) = c1*c1*c1;
273  xi(21) = a; eta(21) = 0.0; zeta(21) = -a; weight(21) = c1*c1*c2;
274  xi(22) = a; eta(22) = 0.0; zeta(22) = 0.0; weight(22) = c1*c2*c2;
275  xi(23) = a; eta(23) = 0.0; zeta(23) = a; weight(23) = c1*c1*c2;
276  xi(24) = a; eta(24) = a; zeta(24) = -a; weight(24) = c1*c1*c1;
277  xi(25) = a; eta(25) = a; zeta(25) = 0.0; weight(25) = c1*c1*c2;
278  xi(26) = a; eta(26) = a; zeta(26) = a; weight(26) = c1*c1*c1;
279 }
280 
281 void gauss_points_tetra1(Epetra_SerialDenseVector & weight, Epetra_SerialDenseVector & xi, Epetra_SerialDenseVector & eta, Epetra_SerialDenseVector & zeta){
282  weight.Resize(1);
283  xi.Resize(1); eta.Resize(1); zeta.Resize(1);
284  weight(0) = 1.0/6.0;
285  xi(0) = 1.0/4.0;
286  eta(0) = 1.0/4.0;
287  zeta(0) = 1.0/4.0;
288 }
289 void gauss_points_tetra4(Epetra_SerialDenseVector & weight, Epetra_SerialDenseVector & xi, Epetra_SerialDenseVector & eta, Epetra_SerialDenseVector & zeta){
290  weight.Resize(4); xi.Resize(4); eta.Resize(4); zeta.Resize(4);
291  double alpha = (5.0 - sqrt(5.0))/20.0;
292  double beta = (5.0 + 3.0*sqrt(5.0))/20.0;
293  weight(0) = 1.0/24.0; weight(1) = 1.0/24.0; weight(2) = 1.0/24.0; weight(3) = 1.0/24.0;
294  xi(0) = alpha; eta(0) = alpha; zeta(0) = alpha;
295  xi(1) = alpha; eta(1) = alpha; zeta(1) = beta;
296  xi(2) = alpha; eta(2) = beta; zeta(2) = alpha;
297  xi(3) = beta; eta(3) = alpha; zeta(3) = alpha;
298 }
299 void gauss_points_tetra5(Epetra_SerialDenseVector & weight, Epetra_SerialDenseVector & xi, Epetra_SerialDenseVector & eta, Epetra_SerialDenseVector & zeta){
300  weight.Resize(5);xi.Resize(5); eta.Resize(5); zeta.Resize(5);
301  xi(0) = 0.25; xi(1) = 0.50; xi(3) = 0.1666666666666667; xi(4) = 0.1666666666666667; xi(5) = 0.1666666666666667;
302  eta(0) = 0.25; eta(1) = 0.1666666666666667; eta(2) = 0.1666666666666667; eta(3) = 0.1666666666666667; eta(4) = 0.50;
303  zeta(0) = 0.25; zeta(1) = 0.1666666666666667; zeta(2) = 0.1666666666666667; zeta(3) = 0.50; zeta(4) = 0.1666666666666667;
304  weight(0) = -0.8000000000000000/6.0; weight(1) = 0.45/6.0; weight(2) = 0.45/6.0; weight(3) = 0.45/6.0; weight(4) =0.45/6.0;
305 
306 }
307 void gauss_points_tetra11(Epetra_SerialDenseVector & weight, Epetra_SerialDenseVector & xi, Epetra_SerialDenseVector & eta, Epetra_SerialDenseVector & zeta){
308  weight.Resize(11); xi.Resize(11); eta.Resize(11); zeta.Resize(11);
309  xi(0) = 0.25; xi(1) = 0.7857142857142857; xi(2) = 0.0714285714285714; xi(3) = 0.0714285714285714; xi(4) = 0.0714285714285714; xi(5) = 0.1005964238332008; xi(6) = 0.3994035761667992; xi(7) = 0.3994035761667992; xi(8) = 0.3994035761667992; xi(9) = 0.1005964238332008; xi(10) = 0.1005964238332008;
310  eta(0) = 0.25; eta(1) = 0.0714285714285714; eta(2) = 0.0714285714285714; eta(3) = 0.0714285714285714; eta(4) = 0.7857142857142857; eta(5) = 0.3994035761667992; eta(6) = 0.1005964238332008; eta(7) = 0.3994035761667992; eta(8) = 0.1005964238332008; eta(9) = 0.3994035761667992; eta(10) = 0.1005964238332008;
311  zeta(0) = 0.25; zeta(1) = 0.0714285714285714; zeta(2) = 0.0714285714285714; zeta(3) = 0.7857142857142857; zeta(4) = 0.0714285714285714; zeta(5) = 0.3994035761667992; zeta(6) = 0.3994035761667992; zeta(7) = 0.1005964238332008; zeta(8) = 0.1005964238332008; zeta(9) = 0.1005964238332008; zeta(10) = 0.3994035761667992;
312  weight(0) = -0.0789333333333333/6.0; weight(1) = 0.0457333333333333/6.0; weight(2) = 0.0457333333333333/6.0; weight(3) = 0.0457333333333333/6.0; weight(4) = 0.0457333333333333/6.0; weight(5) = 0.1493333333333333/6.0; weight(6) = 0.1493333333333333/6.0; weight(7) = 0.1493333333333333/6.0; weight(8) = 0.1493333333333333/6.0; weight(9) = 0.1493333333333333/6.0; weight(10) = 0.1493333333333333/6.0;
313 }
314 
315 int pnpoly(int & nvert, Epetra_SerialDenseVector & vertx, Epetra_SerialDenseVector & verty, double & testx, double & testy){
316  int i, j, c = 0;
317  double a, b;
318  for (i = 0, j = nvert-1; i < nvert; j = i++) {
319  a = (verty[j]-verty[i])/(vertx[j]-vertx[i]);
320  b = verty[i]-a*vertx[i];
321  if (testy - a*testx - b == 0){
322  return 1;
323  }
324  if ( ((verty[i]>testy) != (verty[j]>testy)) &&
325  (testx < (vertx[j]-vertx[i]) * (testy-verty[i]) / (verty[j]-verty[i]) + vertx[i]) ){
326  c = !c;
327  }
328 
329  }
330  return c;
331 }
gauss_points_quad1
void gauss_points_quad1(Epetra_SerialDenseVector &weight, Epetra_SerialDenseVector &xi, Epetra_SerialDenseVector &eta)
Definition: fepp.cpp:202
X
X
Definition: run_meanModels.m:17
j
for j
Definition: costFunction.m:42
quad4::shape_functions
void shape_functions(Epetra_SerialDenseVector &N, double &xi, double &eta)
Definition: fepp.cpp:28
hexa8::shape_functions
void shape_functions(Epetra_SerialDenseVector &N, double &xi, double &eta, double &zeta)
Definition: fepp.cpp:88
gauss_points_quad9
void gauss_points_quad9(Epetra_SerialDenseVector &weight, Epetra_SerialDenseVector &xi, Epetra_SerialDenseVector &eta)
Definition: fepp.cpp:215
quad4::d_shape_functions
void d_shape_functions(Epetra_SerialDenseMatrix &D, double &xi, double &eta)
Definition: fepp.cpp:34
quad4::dX_shape_functions
void dX_shape_functions(Epetra_SerialDenseMatrix &D, Epetra_SerialDenseMatrix &DX, double &jac, Epetra_SerialDenseMatrix &JacobianMatrix)
Definition: fepp.cpp:40
tri6::shape_functions
void shape_functions(Epetra_SerialDenseVector &N, double &xi, double &eta)
Definition: fepp.cpp:51
hexa8::d_shape_functions
void d_shape_functions(Epetra_SerialDenseMatrix &D, double &xi, double &eta, double &zeta)
Definition: fepp.cpp:99
gauss_points_tetra1
void gauss_points_tetra1(Epetra_SerialDenseVector &weight, Epetra_SerialDenseVector &xi, Epetra_SerialDenseVector &eta, Epetra_SerialDenseVector &zeta)
Definition: fepp.cpp:281
dX_shape_functions
void dX_shape_functions(Epetra_SerialDenseMatrix &D, Epetra_SerialDenseMatrix JacobianMatrix, double &jac, Epetra_SerialDenseMatrix &DX)
Definition: fepp.cpp:152
tetra10::d_shape_functions
void d_shape_functions(Epetra_SerialDenseMatrix &D, double &xi, double &eta, double &zeta)
Definition: fepp.cpp:139
jacobian_det
void jacobian_det(Epetra_SerialDenseMatrix &JacobianMatrix, double &jac)
Definition: fepp.cpp:175
gauss_points_tetra5
void gauss_points_tetra5(Epetra_SerialDenseVector &weight, Epetra_SerialDenseVector &xi, Epetra_SerialDenseVector &eta, Epetra_SerialDenseVector &zeta)
Definition: fepp.cpp:299
gauss_points_tri3
void gauss_points_tri3(Epetra_SerialDenseVector &weight, Epetra_SerialDenseVector &xi, Epetra_SerialDenseVector &eta)
Definition: fepp.cpp:189
tri3::dX_shape_functions
void dX_shape_functions(Epetra_SerialDenseMatrix &D, Epetra_SerialDenseMatrix &DX, double &jac, Epetra_SerialDenseMatrix &JacobianMatrix)
Definition: fepp.cpp:17
tetra4::shape_functions
void shape_functions(Epetra_SerialDenseVector &N, double &xi, double &eta, double &zeta)
Definition: fepp.cpp:110
gauss_points_hexa8
void gauss_points_hexa8(Epetra_SerialDenseVector &weight, Epetra_SerialDenseVector &xi, Epetra_SerialDenseVector &eta, Epetra_SerialDenseVector &zeta)
Definition: fepp.cpp:234
gauss_points_tri4
void gauss_points_tri4(Epetra_SerialDenseVector &weight, Epetra_SerialDenseVector &xi, Epetra_SerialDenseVector &eta)
Definition: fepp.cpp:195
jacobian_det_faces
void jacobian_det_faces(Epetra_SerialDenseMatrix &JacobianMatrix, double &jac)
Definition: fepp.cpp:179
tri3::shape_functions
void shape_functions(Epetra_SerialDenseVector &N, double &xi, double &eta)
Definition: fepp.cpp:7
gauss_points_hexa4
void gauss_points_hexa4(Epetra_SerialDenseVector &weight, Epetra_SerialDenseVector &xi, Epetra_SerialDenseVector &eta, Epetra_SerialDenseVector &zeta)
Definition: fepp.cpp:223
gauss_points_tetra4
void gauss_points_tetra4(Epetra_SerialDenseVector &weight, Epetra_SerialDenseVector &xi, Epetra_SerialDenseVector &eta, Epetra_SerialDenseVector &zeta)
Definition: fepp.cpp:289
gauss_points_tri1
void gauss_points_tri1(Epetra_SerialDenseVector &weight, Epetra_SerialDenseVector &xi, Epetra_SerialDenseVector &eta)
Definition: fepp.cpp:183
gauss_points_hexa27
void gauss_points_hexa27(Epetra_SerialDenseVector &weight, Epetra_SerialDenseVector &xi, Epetra_SerialDenseVector &eta, Epetra_SerialDenseVector &zeta)
Definition: fepp.cpp:248
tri6::dX_shape_functions
void dX_shape_functions(Epetra_SerialDenseMatrix &D, Epetra_SerialDenseMatrix &DX, double &jac, Epetra_SerialDenseMatrix &JacobianMatrix)
Definition: fepp.cpp:77
beta
modelParameters beta
Definition: run_station15.m:10
tetra4::d_shape_functions
void d_shape_functions(Epetra_SerialDenseMatrix &D, double &xi, double &eta, double &zeta)
Definition: fepp.cpp:118
tri3::d_shape_functions
void d_shape_functions(Epetra_SerialDenseMatrix &D, double &xi, double &eta)
Definition: fepp.cpp:12
i
for i
Definition: costFunction.m:38
gauss_points_tetra11
void gauss_points_tetra11(Epetra_SerialDenseVector &weight, Epetra_SerialDenseVector &xi, Epetra_SerialDenseVector &eta, Epetra_SerialDenseVector &zeta)
Definition: fepp.cpp:307
gauss_points_quad4
void gauss_points_quad4(Epetra_SerialDenseVector &weight, Epetra_SerialDenseVector &xi, Epetra_SerialDenseVector &eta)
Definition: fepp.cpp:208
jacobian_matrix
void jacobian_matrix(Epetra_SerialDenseMatrix &X, Epetra_SerialDenseMatrix &D, Epetra_SerialDenseMatrix &JacobianMatrix)
Definition: fepp.cpp:171
tri6::d_shape_functions
void d_shape_functions(Epetra_SerialDenseMatrix &D, double &xi, double &eta)
Definition: fepp.cpp:61
pnpoly
int pnpoly(int &nvert, Epetra_SerialDenseVector &vertx, Epetra_SerialDenseVector &verty, double &testx, double &testy)
Definition: fepp.cpp:315
tetra10::shape_functions
void shape_functions(Epetra_SerialDenseVector &N, double &xi, double &eta, double &zeta)
Definition: fepp.cpp:125
eta
output eta
Definition: costFunction.m:33