/[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 616 - (show annotations)
Wed Mar 22 02:46:56 2006 UTC (13 years, 5 months ago) by elspeth
File MIME type: text/plain
File size: 8509 byte(s)
Copyright added to more source files.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26