/[escript]/branches/domexper/dudley/src/ShapeFunctions.c
ViewVC logotype

Contents of /branches/domexper/dudley/src/ShapeFunctions.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3170 - (show annotations)
Thu Sep 9 05:40:25 2010 UTC (8 years, 7 months ago) by jfenwick
File MIME type: text/plain
File size: 7185 byte(s)


1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13
14
15 /**************************************************************/
16
17 /* Dudley: Shape functions */
18
19 /**************************************************************/
20
21 #include "ShapeFunctions.h"
22
23
24
25
26 Dudley_ShapeFunction_Evaluation Dudley_Shape_Point1;
27 Dudley_ShapeFunction_Evaluation Dudley_Shape_Line2;
28 Dudley_ShapeFunction_Evaluation Dudley_Shape_Tri3;
29 Dudley_ShapeFunction_Evaluation Dudley_Shape_Tet4;
30
31
32 Dudley_ShapeFunctionInfo Dudley_ShapeFunction_InfoList[]={
33 {Point1Shape, "Point1", 0, 1, 1, 1, Dudley_Shape_Point1 } ,
34 {Line2Shape, "Line2", 1, 2, 1, 2, Dudley_Shape_Line2 } ,
35 {Tri3Shape, "Tri3", 2, 3, 1, 3, Dudley_Shape_Tri3 },
36 {Tet4Shape, "Tet4", 3, 4, 1, 4, Dudley_Shape_Tet4 },
37 {NoShape, "NoType", 0, 1, 1, 1, Dudley_Shape_Point1 }
38 };
39
40
41 /******************************************************************************************************************************
42
43 creates an evaluation of the ShapeFunction on the given quadrature scheme.
44 if the spatial dimension of the scheme and the shape functions don't match
45
46 if QuadNodes==Null or QuadWeights==Null the shape functions method is used to generate a quadrature scheme with numQuasNodes
47 nodes. otherwise its assumed that a quadraure scheme is given on these array and copy is created within the structure.
48
49 */
50 Dudley_ShapeFunction* Dudley_ShapeFunction_alloc(Dudley_ShapeFunctionTypeId id,int numQuadDim, int numQuadNodes, double *QuadNodes, double *QuadWeights) {
51 Dudley_ShapeFunction *out=NULL;
52 int numDim, numShapes, i, q;
53
54 numDim=Dudley_ShapeFunction_InfoList[id].numDim;
55 numShapes=Dudley_ShapeFunction_InfoList[id].numShapes;
56
57 if (numQuadDim>numDim) {
58 Dudley_setError(VALUE_ERROR,"Dudley_ShapeFunction_alloc: spatial dimension of quadrature scheme is bigger then spatial dimension of shape function.");
59 return NULL;
60 }
61
62 /* allocate the Dudley_ShapeFunction to be returned: */
63
64 out=MEMALLOC(1,Dudley_ShapeFunction);
65 if (Dudley_checkPtr(out)) return NULL;
66
67
68 out->Type=Dudley_ShapeFunction_getInfo(id);
69 out->numQuadNodes=numQuadNodes;
70 out->QuadNodes=NULL;
71 out->QuadWeights=NULL;
72 out->S=NULL;
73 out->dSdv=NULL;
74 out->reference_counter=0;
75
76 /* allocate memory: */
77
78 out->QuadNodes=MEMALLOC(numQuadNodes*numDim,double);
79 out->QuadWeights=MEMALLOC(numQuadNodes,double);
80 out->S=MEMALLOC(numShapes*numQuadNodes,double);
81 out->dSdv=MEMALLOC(numShapes*numDim*numQuadNodes,double);
82 if ( Dudley_checkPtr(out->QuadNodes) || Dudley_checkPtr(out->QuadWeights) || Dudley_checkPtr(out->S) || Dudley_checkPtr(out->dSdv) ) {
83 Dudley_ShapeFunction_dealloc(out);
84 return NULL;
85 }
86
87 /* set the quadrature nodes (missing values are filled with 0): */
88
89 for (q=0;q<numQuadNodes;q++) {
90 for (i=0;i<numQuadDim;i++) out->QuadNodes[INDEX2(i,q,numDim)]=QuadNodes[INDEX2(i,q,numQuadDim)];
91 for (i=numQuadDim;i<numDim;i++) out->QuadNodes[INDEX2(i,q,numDim)]=0;
92 out->QuadWeights[q]=QuadWeights[q];
93 }
94
95 /* eval shape functions on quadrature node: */
96
97 out->Type->getValues(numQuadNodes,out->QuadNodes,out->S,out->dSdv);
98
99 if (! Dudley_noError()) {
100 Dudley_ShapeFunction_dealloc(out);
101 return NULL;
102 }
103
104 /* all done: */
105 out->reference_counter=1;
106 return out;
107 }
108
109 Dudley_ShapeFunction* Dudley_ShapeFunction_reference(Dudley_ShapeFunction* in) {
110 if (in!=NULL) ++(in->reference_counter);
111 return in;
112 }
113 /**************************************************************/
114
115 void Dudley_ShapeFunction_dealloc(Dudley_ShapeFunction* in) {
116 if (in!=NULL) {
117 in->reference_counter--;
118 if (in->reference_counter<1) {
119 MEMFREE(in->QuadNodes);
120 MEMFREE(in->QuadWeights);
121 MEMFREE(in->S);
122 MEMFREE(in->dSdv);
123 MEMFREE(in);
124 }
125 }
126 }
127
128 /**************************************************************/
129
130 Dudley_ShapeFunctionTypeId Dudley_ShapeFunction_getTypeId(char* element_type)
131 {
132 int ptr=0;
133 Dudley_ShapeFunctionTypeId out=NoShape;
134 while (Dudley_ShapeFunction_InfoList[ptr].TypeId!=NoShape && out==NoShape) {
135 if (strcmp(element_type,Dudley_ShapeFunction_InfoList[ptr].Name)==0) out=Dudley_ShapeFunction_InfoList[ptr].TypeId;
136 ptr++;
137 }
138 return out;
139 }
140
141 Dudley_ShapeFunctionInfo* Dudley_ShapeFunction_getInfo(Dudley_ShapeFunctionTypeId id)
142 {
143 int ptr=0;
144 Dudley_ShapeFunctionInfo* out=NULL;
145 while (Dudley_ShapeFunction_InfoList[ptr].TypeId!=NoShape && out==NULL) {
146 if (Dudley_ShapeFunction_InfoList[ptr].TypeId==id) out=&(Dudley_ShapeFunction_InfoList[ptr]);
147 ptr++;
148 }
149 if (out==NULL) {
150 Dudley_setError(VALUE_ERROR,"Dudley_ShapeFunctionInfo_getInfo: canot find requested shape function");
151 }
152 return out;
153 }
154
155 /**************************************************************/
156
157 #define V(_K_,_I_) v[INDEX2((_K_)-1,(_I_),DIM)]
158 #define S(_J_,_I_) s[S_INDEX((_J_)-1,(_I_),NUMSHAPES)]
159 #define DSDV(_J_,_K_,_I_) dsdv[DSDV_INDEX((_J_)-1,(_K_)-1,(_I_),NUMSHAPES,DIM)]
160
161 /**************************************************************/
162
163 void Dudley_Shape_Point1(int NumV,double* v,double* s,double* dsdv) {
164 #define NUMSHAPES 1
165 #define DIM 0
166 int i;
167 #pragma ivdep
168 for (i=0;i<NumV;i++) {
169 S(1,i)=1.;
170 }
171 #undef NUMSHAPES
172 #undef DIM
173 }
174
175 /**************************************************************/
176
177 void Dudley_Shape_Line2(int NumV,double* v,double* s,double* dsdv) {
178 #define NUMSHAPES 2
179 #define DIM 1
180 register double x;
181 int i;
182 #pragma ivdep
183 for (i=0;i<NumV;i++) {
184 x=V(1,i);
185 S(1,i)=1.-x;
186 S(2,i)= x;
187 DSDV(1,1,i)=-1.;
188 DSDV(2,1,i)= 1.;
189 }
190 #undef NUMSHAPES
191 #undef DIM
192 }
193
194 /**************************************************************/
195
196 void Dudley_Shape_Tri3(int NumV,double* v,double* s,double* dsdv) {
197 #define NUMSHAPES 3
198 #define DIM 2
199 register double x,y;
200 int i;
201 #pragma ivdep
202 for (i=0;i<NumV;i++) {
203 x=V(1,i);
204 y=V(2,i);
205 S(1,i)=1.-x-y;
206 S(2,i)= x;
207 S(3,i)= y;
208 DSDV(1,1,i)=-1.;
209 DSDV(1,2,i)=-1.;
210 DSDV(2,1,i)= 1.;
211 DSDV(2,2,i)= 0.;
212 DSDV(3,1,i)= 0.;
213 DSDV(3,2,i)= 1.;
214 }
215 #undef NUMSHAPES
216 #undef DIM
217 }
218
219 /**************************************************************/
220
221 void Dudley_Shape_Tet4(int NumV,double* v,double* s,double* dsdv) {
222 #define NUMSHAPES 4
223 #define DIM 3
224 register double x,y,z;
225 int i;
226 #pragma ivdep
227 for (i=0;i<NumV;i++) {
228 x=V(1,i);
229 y=V(2,i);
230 z=V(3,i);
231 S(1,i)=1.-x-y-z;
232 S(2,i)=x;
233 S(3,i)=y;
234 S(4,i)=z;
235 DSDV(1,1,i)=-1.;
236 DSDV(1,2,i)=-1.;
237 DSDV(1,3,i)=-1.;
238 DSDV(2,1,i)= 1.;
239 DSDV(2,2,i)= 0.;
240 DSDV(2,3,i)= 0.;
241 DSDV(3,1,i)= 0.;
242 DSDV(3,2,i)= 1.;
243 DSDV(3,3,i)= 0.;
244 DSDV(4,1,i)= 0.;
245 DSDV(4,2,i)= 0.;
246 DSDV(4,3,i)= 1.;
247 }
248 #undef NUMSHAPES
249 #undef DIM
250 }
251
252 #undef V
253 #undef S
254 #undef DSDV
255
256

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26