/[escript]/trunk/dudley/src/IndexList.cpp
ViewVC logotype

Annotation of /trunk/dudley/src/IndexList.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2748 - (hide annotations)
Tue Nov 17 07:32:59 2009 UTC (9 years, 10 months ago) by gross
Original Path: trunk/finley/src/IndexList.c
File MIME type: text/plain
File size: 8838 byte(s)
Macro elements are implemented now. VTK writer for macro elements still needs testing.
1 jgs 82
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4 jfenwick 2548 * Copyright (c) 2003-2009 by University of Queensland
5 ksteube 1811 * 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 jgs 82
14 ksteube 1811
15 jgs 82 /**************************************************************/
16    
17 ksteube 1312 /* Finley: Converting an element list into a matrix shape */
18 jgs 82
19     /**************************************************************/
20    
21     #include "IndexList.h"
22    
23 ksteube 1312 /* Translate from distributed/local array indices to global indices */
24    
25 jgs 82 /**************************************************************/
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 ksteube 971 void Finley_IndexList_insertElements(Finley_IndexList* index_list, Finley_ElementFile* elements,
31 ksteube 1312 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 gross 2748 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 jgs 82 if (elements!=NULL) {
38 ksteube 1312 NN=elements->numNodes;
39 gross 2748 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 ksteube 1312 for (e=0;e<elements->numElements;e++) {
64     if (elements->Color[e]==color) {
65 gross 2748 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 ksteube 1312 }
75     }
76     }
77 jgs 82 }
78     return;
79     }
80    
81 gross 2748
82 ksteube 1312 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 gross 2748 /* this does not resolve macro elements */
86     index_t color;
87 phornby 1628 dim_t e,kr,kc,icol,irow, NN;
88 ksteube 1312 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 gross 1722 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 gross 2748 /* this does not resolve macro elements */
113 gross 1722 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 ksteube 1312
137 jgs 82 /* inserts row index row into the Finley_IndexList in if it does not exist */
138    
139 jgs 123 void Finley_IndexList_insertIndex(Finley_IndexList* in, index_t index) {
140     dim_t i;
141 jgs 82 /* is index in in? */
142     for (i=0;i<in->n;i++) {
143 jgs 102 if (in->index[i]==index) return;
144 jgs 82 }
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 jgs 102 in->extension=TMPMEMALLOC(1,Finley_IndexList);
150 jgs 82 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 ksteube 1312 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 jgs 82 if (in==NULL) {
169     return 0;
170     } else {
171 ksteube 1312 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 jgs 82 }
177     }
178    
179     /* count the number of row indices in the Finley_IndexList in */
180    
181 ksteube 1312 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 jgs 82 if (in!=NULL) {
185 ksteube 1312 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 jgs 82 }
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 ksteube 1312 /* creates a Paso_pattern from a range of indices */
208 gross 1552 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 ksteube 1312 {
210     dim_t *ptr=NULL;
211     register dim_t s,i,itmp;
212     index_t *index=NULL;
213     Paso_Pattern* out=NULL;
214    
215 gross 1552 ptr=MEMALLOC(n+1-n0,index_t);
216 ksteube 1312 if (! Finley_checkPtr(ptr) ) {
217     /* get the number of connections per row */
218     #pragma omp parallel for schedule(static) private(i)
219 gross 1552 for(i=n0;i<n;++i) {
220     ptr[i-n0]=Finley_IndexList_count(&index_list[i],range_min,range_max);
221 ksteube 1312 }
222     /* accumulate ptr */
223     s=0;
224 gross 1552 for(i=n0;i<n;++i) {
225     itmp=ptr[i-n0];
226     ptr[i-n0]=s;
227 ksteube 1312 s+=itmp;
228     }
229 gross 1552 ptr[n-n0]=s;
230 ksteube 1312 /* fill index */
231 gross 1552 index=MEMALLOC(ptr[n-n0],index_t);
232 ksteube 1312 if (! Finley_checkPtr(index)) {
233     #pragma omp parallel for schedule(static)
234 gross 1552 for(i=n0;i<n;++i) {
235     Finley_IndexList_toArray(&index_list[i],&index[ptr[i-n0]],range_min,range_max,index_offset);
236 ksteube 1312 }
237 gross 2551 out=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,n-n0,range_max+index_offset,ptr,index);
238 ksteube 1312 }
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