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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2748 - (show annotations)
Tue Nov 17 07:32:59 2009 UTC (9 years, 7 months ago) by gross
File MIME type: text/plain
File size: 8838 byte(s)
Macro elements are implemented now. VTK writer for macro elements still needs testing.
1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2009 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 /* Finley: Converting an element list into a matrix shape */
18
19 /**************************************************************/
20
21 #include "IndexList.h"
22
23 /* Translate from distributed/local array indices to global indices */
24
25 /**************************************************************/
26 /* inserts the contributions from the element matrices of elements
27 into the row index col. If symmetric is set, only the upper
28 triangle of the matrix is stored. */
29
30 void Finley_IndexList_insertElements(Finley_IndexList* index_list, Finley_ElementFile* elements,
31 bool_t reduce_row_order, index_t* row_map,
32 bool_t reduce_col_order, index_t* col_map) {
33 /* index_list is an array of linked lists. Each entry is a row (DOF) and contains the indices to the non-zero columns */
34 index_t color;
35 Finley_ReferenceElement*refElement;
36 dim_t e, kr, kc, NN_row, NN_col, icol, irow, NN, *row_node=NULL, *col_node=NULL, isub, numSub;
37 if (elements!=NULL) {
38 NN=elements->numNodes;
39 refElement= Finley_ReferenceElementSet_borrowReferenceElement(elements->referenceElementSet, FALSE);
40
41 if (reduce_col_order) {
42 numSub=1;
43 col_node=refElement->Type->linearNodes;
44 NN_col=(refElement->LinearBasisFunctions->Type->numShapes) * (refElement->Type->numSides) ;
45 } else {
46 numSub=refElement->Type->numSubElements;
47 col_node=refElement->Type->subElementNodes;
48 NN_col=(refElement->BasisFunctions->Type->numShapes) * (refElement->Type->numSides) ;
49 }
50
51 if (reduce_row_order) {
52 numSub=1;
53 row_node=refElement->Type->linearNodes;
54 NN_row=(refElement->LinearBasisFunctions->Type->numShapes) * (refElement->Type->numSides) ;
55 } else {
56 numSub=refElement->Type->numSubElements;
57 row_node=refElement->Type->subElementNodes;
58 NN_row=(refElement->BasisFunctions->Type->numShapes) * (refElement->Type->numSides) ;
59 }
60
61 for (color=elements->minColor;color<=elements->maxColor;color++) {
62 #pragma omp for private(e,irow,kr,kc,icol,isub) schedule(static)
63 for (e=0;e<elements->numElements;e++) {
64 if (elements->Color[e]==color) {
65 for (isub=0;isub<numSub; isub++) {
66 for (kr=0;kr<NN_row;kr++) {
67 irow=row_map[elements->Nodes[INDEX2(row_node[INDEX2(kr,isub,NN_row)],e,NN)]];
68 for (kc=0;kc<NN_col;kc++) {
69 icol=col_map[elements->Nodes[INDEX2(col_node[INDEX2(kc,isub,NN_col)],e,NN)]];
70 Finley_IndexList_insertIndex(&(index_list[irow]),icol);
71 }
72 }
73 }
74 }
75 }
76 }
77 }
78 return;
79 }
80
81
82 void Finley_IndexList_insertElementsWithRowRange(Finley_IndexList* index_list, index_t firstRow, index_t lastRow,
83 Finley_ElementFile* elements, index_t* row_map, index_t* col_map)
84 {
85 /* this does not resolve macro elements */
86 index_t color;
87 dim_t e,kr,kc,icol,irow, NN;
88 if (elements!=NULL) {
89 NN=elements->numNodes;
90 for (color=elements->minColor;color<=elements->maxColor;color++) {
91 #pragma omp for private(e,irow,kr,kc,icol) schedule(static)
92 for (e=0;e<elements->numElements;e++) {
93 if (elements->Color[e]==color) {
94 for (kr=0;kr<NN;kr++) {
95 irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];
96 if ((firstRow<=irow) && (irow < lastRow)) {
97 irow-=firstRow;
98 for (kc=0;kc<NN;kc++) {
99 icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];
100 Finley_IndexList_insertIndex(&(index_list[irow]),icol);
101 }
102 }
103 }
104 }
105 }
106 }
107 }
108 }
109 void Finley_IndexList_insertElementsWithRowRangeNoMainDiagonal(Finley_IndexList* index_list, index_t firstRow, index_t lastRow,
110 Finley_ElementFile* elements, index_t* row_map, index_t* col_map)
111 {
112 /* this does not resolve macro elements */
113 index_t color;
114 dim_t e,kr,kc,icol,irow, NN,irow_loc;
115 if (elements!=NULL) {
116 NN=elements->numNodes;
117 for (color=elements->minColor;color<=elements->maxColor;color++) {
118 #pragma omp for private(e,irow,kr,kc,icol,irow_loc) schedule(static)
119 for (e=0;e<elements->numElements;e++) {
120 if (elements->Color[e]==color) {
121 for (kr=0;kr<NN;kr++) {
122 irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];
123 if ((firstRow<=irow) && (irow < lastRow)) {
124 irow_loc=irow-firstRow;
125 for (kc=0;kc<NN;kc++) {
126 icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];
127 if (icol != irow) Finley_IndexList_insertIndex(&(index_list[irow_loc]),icol);
128 }
129 }
130 }
131 }
132 }
133 }
134 }
135 }
136
137 /* inserts row index row into the Finley_IndexList in if it does not exist */
138
139 void Finley_IndexList_insertIndex(Finley_IndexList* in, index_t index) {
140 dim_t i;
141 /* is index in in? */
142 for (i=0;i<in->n;i++) {
143 if (in->index[i]==index) return;
144 }
145 /* index could not be found */
146 if (in->n==INDEXLIST_LENGTH) {
147 /* if in->index is full check the extension */
148 if (in->extension==NULL) {
149 in->extension=TMPMEMALLOC(1,Finley_IndexList);
150 if (Finley_checkPtr(in->extension)) return;
151 in->extension->n=0;
152 in->extension->extension=NULL;
153 }
154 Finley_IndexList_insertIndex(in->extension,index);
155 } else {
156 /* insert index into in->index*/
157 in->index[in->n]=index;
158 in->n++;
159 }
160 }
161
162 /* counts the number of row indices in the Finley_IndexList in */
163
164 dim_t Finley_IndexList_count(Finley_IndexList* in, index_t range_min,index_t range_max) {
165 dim_t i;
166 dim_t out=0;
167 register index_t itmp;
168 if (in==NULL) {
169 return 0;
170 } else {
171 for (i=0;i<in->n;i++) {
172 itmp=in->index[i];
173 if ((itmp>=range_min) && (range_max>itmp)) ++out;
174 }
175 return out+Finley_IndexList_count(in->extension, range_min,range_max);
176 }
177 }
178
179 /* count the number of row indices in the Finley_IndexList in */
180
181 void Finley_IndexList_toArray(Finley_IndexList* in, index_t* array, index_t range_min,index_t range_max, index_t index_offset) {
182 dim_t i, ptr;
183 register index_t itmp;
184 if (in!=NULL) {
185 ptr=0;
186 for (i=0;i<in->n;i++) {
187 itmp=in->index[i];
188 if ((itmp>=range_min) && (range_max>itmp)) {
189 array[ptr]=itmp+index_offset;
190 ptr++;
191 }
192
193 }
194 Finley_IndexList_toArray(in->extension,&(array[ptr]), range_min, range_max, index_offset);
195 }
196 }
197
198 /* deallocates the Finley_IndexList in by recursive calls */
199
200 void Finley_IndexList_free(Finley_IndexList* in) {
201 if (in!=NULL) {
202 Finley_IndexList_free(in->extension);
203 TMPMEMFREE(in);
204 }
205 }
206
207 /* creates a Paso_pattern from a range of indices */
208 Paso_Pattern* Finley_IndexList_createPattern(dim_t n0, dim_t n,Finley_IndexList* index_list,index_t range_min,index_t range_max,index_t index_offset)
209 {
210 dim_t *ptr=NULL;
211 register dim_t s,i,itmp;
212 index_t *index=NULL;
213 Paso_Pattern* out=NULL;
214
215 ptr=MEMALLOC(n+1-n0,index_t);
216 if (! Finley_checkPtr(ptr) ) {
217 /* get the number of connections per row */
218 #pragma omp parallel for schedule(static) private(i)
219 for(i=n0;i<n;++i) {
220 ptr[i-n0]=Finley_IndexList_count(&index_list[i],range_min,range_max);
221 }
222 /* accumulate ptr */
223 s=0;
224 for(i=n0;i<n;++i) {
225 itmp=ptr[i-n0];
226 ptr[i-n0]=s;
227 s+=itmp;
228 }
229 ptr[n-n0]=s;
230 /* fill index */
231 index=MEMALLOC(ptr[n-n0],index_t);
232 if (! Finley_checkPtr(index)) {
233 #pragma omp parallel for schedule(static)
234 for(i=n0;i<n;++i) {
235 Finley_IndexList_toArray(&index_list[i],&index[ptr[i-n0]],range_min,range_max,index_offset);
236 }
237 out=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,n-n0,range_max+index_offset,ptr,index);
238 }
239 }
240 if (! Finley_noError()) {
241 MEMFREE(ptr);
242 MEMFREE(index);
243 Paso_Pattern_free(out);
244 }
245 return out;
246 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26