/[escript]/trunk/finley/src/Mesh_rec4.c
ViewVC logotype

Contents of /trunk/finley/src/Mesh_rec4.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 149 - (show annotations)
Thu Sep 1 03:31:39 2005 UTC (14 years, 2 months ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Mesh_rec4.c
File MIME type: text/plain
File size: 8441 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-09-01

1 /**************************************************************/
2
3 /* Finley: generates rectangular meshes */
4
5 /* Generates a numElements[0] x numElements[1] mesh with first order elements (Rec4) in the rectangle */
6 /* [0,Length[0]] x [0,Length[1]]. order is the desired accuracy of the integration scheme. */
7
8
9 /**************************************************************/
10
11 /* Copyrights by ACcESS Australia 2003 */
12 /* Author: gross@access.edu.au */
13 /* Version: $Id$ */
14
15 /**************************************************************/
16
17 #include "Common.h"
18 #include "Finley.h"
19 #include "Mesh.h"
20 #include "RectangularMesh.h"
21
22 /**************************************************************/
23
24 Finley_Mesh* Finley_RectangularMesh_Rec4(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace) {
25 dim_t N0,N1,NE0,NE1,i0,i1,totalNECount,faceNECount,NDOF0,NDOF1,NFaceElements,NUMNODES;
26 index_t k,node0;
27 Finley_Mesh* out;
28 char name[50];
29 double time0=Finley_timer();
30 NE0=MAX(1,numElements[0]);
31 NE1=MAX(1,numElements[1]);
32 N0=NE0+1;
33 N1=NE1+1;
34
35 NFaceElements=0;
36 if (!periodic[0]) {
37 NDOF0=N0;
38 NFaceElements+=2*NE1;
39 } else {
40 NDOF0=N0-1;
41 }
42 if (!periodic[1]) {
43 NDOF1=N1;
44 NFaceElements+=2*NE0;
45 } else {
46 NDOF1=N1-1;
47 }
48
49 /* allocate mesh: */
50
51 sprintf(name,"Rectangular %d x %d mesh",N0,N1);
52 out=Finley_Mesh_alloc(name,2,order);
53 if (Finley_ErrorCode!=NO_ERROR) return NULL;
54
55 out->Elements=Finley_ElementFile_alloc(Rec4,out->order);
56 if (useElementsOnFace) {
57 out->FaceElements=Finley_ElementFile_alloc(Rec4Face,out->order);
58 out->ContactElements=Finley_ElementFile_alloc(Rec4Face_Contact,out->order);
59 } else {
60 out->FaceElements=Finley_ElementFile_alloc(Line2,out->order);
61 out->ContactElements=Finley_ElementFile_alloc(Line2_Contact,out->order);
62 }
63 out->Points=Finley_ElementFile_alloc(Point1,out->order);
64 if (Finley_ErrorCode!=NO_ERROR) {
65 Finley_Mesh_dealloc(out);
66 return NULL;
67 }
68
69 /* allocate tables: */
70
71 Finley_NodeFile_allocTable(out->Nodes,N0*N1);
72 Finley_ElementFile_allocTable(out->Elements,NE0*NE1);
73 Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
74 if (Finley_ErrorCode!=NO_ERROR) {
75 Finley_Mesh_dealloc(out);
76 return NULL;
77 }
78
79 /* set nodes: */
80
81 #pragma omp parallel for private(i0,i1,k)
82 for (i1=0;i1<N1;i1++) {
83 for (i0=0;i0<N0;i0++) {
84 k=i0+N0*i1;
85 out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE(i0)/DBLE(N0-1)*Length[0];
86 out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];
87 out->Nodes->Id[k]=k;
88 out->Nodes->Tag[k]=0;
89 out->Nodes->degreeOfFreedom[k]=(i0%NDOF0) +N0*(i1%NDOF1);
90 }
91 }
92 /* tags for the faces: */
93 if (!periodic[1]) {
94 for (i0=0;i0<N0;i0++) {
95 out->Nodes->Tag[i0+N0*0]+=10;
96 out->Nodes->Tag[i0+N0*(N1-1)]+=20;
97 }
98 }
99 if (!periodic[0]) {
100 for (i1=0;i1<N1;i1++) {
101 out->Nodes->Tag[0+N0*i1]+=1;
102 out->Nodes->Tag[(N0-1)+N0*i1]+=2;
103 }
104 }
105
106 /* set the elements: */
107
108 #pragma omp parallel for private(i0,i1,k,node0)
109 for (i1=0;i1<NE1;i1++) {
110 for (i0=0;i0<NE0;i0++) {
111 k=i0+NE0*i1;
112 node0=i0+i1*N0;
113
114 out->Elements->Id[k]=k;
115 out->Elements->Tag[k]=0;
116 out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1);
117
118 out->Elements->Nodes[INDEX2(0,k,4)]=node0;
119 out->Elements->Nodes[INDEX2(1,k,4)]=node0+1;
120 out->Elements->Nodes[INDEX2(2,k,4)]=node0+N0+1;
121 out->Elements->Nodes[INDEX2(3,k,4)]=node0+N0;
122
123 }
124 }
125 out->Elements->minColor=0;
126 out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0);
127
128 /* face elements: */
129
130 if (useElementsOnFace) {
131 NUMNODES=4;
132 } else {
133 NUMNODES=2;
134 }
135 totalNECount=NE0*NE1;
136 faceNECount=0;
137
138 /* elements on boundary 010 (x2=0): */
139 if (!periodic[1]) {
140 #pragma omp parallel for private(i0,k,node0)
141 for (i0=0;i0<NE0;i0++) {
142 k=i0+faceNECount;
143 node0=i0;
144
145 out->FaceElements->Id[k]=i0+totalNECount;
146 out->FaceElements->Tag[k]=10;
147 out->FaceElements->Color[k]=i0%2;
148
149 if (useElementsOnFace) {
150 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
151 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
152 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;
153 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;
154 } else {
155 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
156 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
157 }
158 }
159 totalNECount+=NE0;
160 faceNECount+=NE0;
161
162 /* ** elements on boundary 020 (x2=1): */
163
164 #pragma omp parallel for private(i0,k,node0)
165 for (i0=0;i0<NE0;i0++) {
166 k=i0+faceNECount;
167 node0=i0+(NE1-1)*N0;
168
169 out->FaceElements->Id[k]=i0+totalNECount;
170 out->FaceElements->Tag[k]=20;
171 out->FaceElements->Color[k]=i0%2+2;
172
173 if (useElementsOnFace) {
174 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0+1;
175 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
176 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0;
177 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
178 } else {
179 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0+1;
180 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
181 }
182 }
183 totalNECount+=NE0;
184 faceNECount+=NE0;
185 }
186 if (!periodic[0]) {
187 /* ** elements on boundary 001 (x1=0): */
188 #pragma omp parallel for private(i1,k,node0)
189 for (i1=0;i1<NE1;i1++) {
190 k=i1+faceNECount;
191 node0=i1*N0;
192
193 out->FaceElements->Id[k]=i1+totalNECount;
194 out->FaceElements->Tag[k]=1;
195 out->FaceElements->Color[k]=i1%2+4;
196
197 if (useElementsOnFace) {
198 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
199 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
200 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;
201 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;
202 } else {
203 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
204 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
205
206 }
207 }
208 totalNECount+=NE1;
209 faceNECount+=NE1;
210
211 /* ** elements on boundary 002 (x1=1): */
212
213 #pragma omp parallel for private(i1,k,node0)
214 for (i1=0;i1<NE1;i1++) {
215 k=i1+faceNECount;
216 node0=(NE0-1)+i1*N0;
217
218 out->FaceElements->Id[k]=i1+totalNECount;
219 out->FaceElements->Tag[k]=2;
220 out->FaceElements->Color[k]=i1%2+6;
221
222 if (useElementsOnFace) {
223 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
224 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
225 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0;
226 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
227 } else {
228 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
229 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
230 }
231 }
232 totalNECount+=NE1;
233 faceNECount+=NE1;
234 }
235 out->FaceElements->minColor=0;
236 out->FaceElements->maxColor=7;
237
238 /* face elements done: */
239
240 /* condense the nodes: */
241
242 Finley_Mesh_resolveNodeIds(out);
243
244 /* prepare mesh for further calculatuions:*/
245
246 Finley_Mesh_prepare(out) ;
247
248 #ifdef Finley_TRACE
249 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
250 #endif
251
252 if (Finley_ErrorCode!=NO_ERROR) {
253 Finley_Mesh_dealloc(out);
254 return NULL;
255 }
256 return out;
257 }
258
259 /*
260 * $Log$
261 * Revision 1.3 2005/09/01 03:31:36 jgs
262 * Merge of development branch dev-02 back to main trunk on 2005-09-01
263 *
264 * Revision 1.2.2.1 2005/08/24 02:02:18 gross
265 * timing output switched off. solver output can be swiched through getSolution(verbose=True) now.
266 *
267 * Revision 1.2 2005/07/08 04:07:54 jgs
268 * Merge of development branch back to main trunk on 2005-07-08
269 *
270 * Revision 1.1.1.1.2.1 2005/06/29 02:34:53 gross
271 * some changes towards 64 integers in finley
272 *
273 * Revision 1.1.1.1 2004/10/26 06:53:57 jgs
274 * initial import of project esys2
275 *
276 * Revision 1.1.1.1 2004/06/24 04:00:40 johng
277 * Initial version of eys using boost-python.
278 *
279 *
280 */
281

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26