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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4496 - (show annotations)
Mon Jul 15 06:53:44 2013 UTC (6 years ago) by caltinay
File size: 6558 byte(s)
finley (WIP):
-moved all of finley into its namespace
-introduced some shared pointers
-Mesh is now a class
-other bits and pieces...

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2013 by University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development since 2012 by School of Earth Sciences
13 *
14 *****************************************************************************/
15
16
17 /****************************************************************************
18
19 Finley: Converting an element list into a matrix shape
20
21 *****************************************************************************/
22
23 #include "IndexList.h"
24 #include "ElementFile.h"
25
26 #include "paso/SystemMatrixPattern.h"
27
28 /* Translate from distributed/local array indices to global indices */
29
30 /****************************************************************************/
31 /* inserts the contributions from the element matrices of elements
32 into the row index col. If symmetric is set, only the upper
33 triangle of the matrix is stored.
34 */
35
36 namespace finley {
37
38 void IndexList_insertElements(IndexList* index_list, ElementFile* elements,
39 bool reduce_row_order, const int* row_map,
40 bool reduce_col_order, const int* col_map)
41 {
42 // index_list is an array of linked lists. Each entry is a row (DOF) and
43 // contains the indices to the non-zero columns
44 if (!elements)
45 return;
46
47 const int NN=elements->numNodes;
48 const_ReferenceElement_ptr refElement(elements->referenceElementSet->
49 borrowReferenceElement(false));
50
51 int NN_row, NN_col, numSub;
52 const int *row_node=NULL, *col_node=NULL;
53 if (reduce_col_order) {
54 numSub=1;
55 col_node=refElement->Type->linearNodes;
56 NN_col=refElement->LinearBasisFunctions->Type->numShapes * refElement->Type->numSides;
57 } else {
58 numSub=refElement->Type->numSubElements;
59 col_node=refElement->Type->subElementNodes;
60 NN_col=refElement->BasisFunctions->Type->numShapes * refElement->Type->numSides;
61 }
62
63 if (reduce_row_order) {
64 numSub=1;
65 row_node=refElement->Type->linearNodes;
66 NN_row=refElement->LinearBasisFunctions->Type->numShapes * refElement->Type->numSides;
67 } else {
68 numSub=refElement->Type->numSubElements;
69 row_node=refElement->Type->subElementNodes;
70 NN_row=refElement->BasisFunctions->Type->numShapes * refElement->Type->numSides;
71 }
72
73 for (int color=elements->minColor; color<=elements->maxColor; color++) {
74 #pragma omp for
75 for (int e=0; e<elements->numElements; e++) {
76 if (elements->Color[e]==color) {
77 for (int isub=0; isub<numSub; isub++) {
78 for (int kr=0; kr<NN_row; kr++) {
79 const int irow=row_map[elements->Nodes[INDEX2(row_node[INDEX2(kr,isub,NN_row)],e,NN)]];
80 for (int kc=0; kc<NN_col; kc++) {
81 const int icol=col_map[elements->Nodes[INDEX2(col_node[INDEX2(kc,isub,NN_col)],e,NN)]];
82 IndexList_insertIndex(index_list[irow], icol);
83 }
84 }
85 }
86 }
87 }
88 }
89 }
90
91 void IndexList_insertElementsWithRowRangeNoMainDiagonal(IndexList* index_list,
92 int firstRow, int lastRow, ElementFile* elements,
93 int* row_map, int* col_map)
94 {
95 if (!elements)
96 return;
97
98 // this does not resolve macro elements
99 const int NN=elements->numNodes;
100 for (int color=elements->minColor; color<=elements->maxColor; color++) {
101 #pragma omp for
102 for (int e=0; e<elements->numElements; e++) {
103 if (elements->Color[e]==color) {
104 for (int kr=0; kr<NN; kr++) {
105 const int irow=row_map[elements->Nodes[INDEX2(kr,e,NN)]];
106 if (firstRow<=irow && irow<lastRow) {
107 const int irow_loc=irow-firstRow;
108 for (int kc=0; kc<NN; kc++) {
109 const int icol=col_map[elements->Nodes[INDEX2(kc,e,NN)]];
110 if (icol != irow)
111 IndexList_insertIndex(index_list[irow_loc], icol);
112 }
113 }
114 }
115 }
116 }
117 }
118 }
119
120 /// inserts a row index into the IndexList in if it does not exist
121 void IndexList_insertIndex(IndexList& in, int index)
122 {
123 if (std::find(in.begin(), in.end(), index) == in.end())
124 in.push_back(index);
125 }
126
127 /// counts the number of row indices in the IndexList in
128 int IndexList_count(const IndexList& in, int range_min, int range_max)
129 {
130 int out=0;
131 for (IndexList::const_iterator it=in.begin(); it!=in.end(); it++) {
132 if (*it >= range_min && range_max > *it)
133 ++out;
134 }
135 return out;
136 }
137
138 /// count the number of row indices in the IndexList in
139 void IndexList_toArray(const IndexList& in, int* array, int range_min, int range_max, int index_offset)
140 {
141 int idx=0;
142 for (IndexList::const_iterator it=in.begin(); it!=in.end(); it++) {
143 if (*it >= range_min && range_max > *it) {
144 array[idx]=(*it)+index_offset;
145 ++idx;
146 }
147 }
148 }
149
150 /// creates a Paso_pattern from a range of indices
151 Paso_Pattern* IndexList_createPattern(int n0, int n,
152 const IndexList* index_list,
153 int range_min, int range_max,
154 int index_offset)
155 {
156 int *ptr=new int[n+1-n0];
157 // get the number of connections per row
158 #pragma omp parallel for
159 for (int i=n0; i<n; ++i) {
160 ptr[i-n0]=IndexList_count(index_list[i], range_min, range_max);
161 }
162 // accumulate ptr
163 int s=0;
164 for (int i=n0; i<n; ++i) {
165 const int itmp=ptr[i-n0];
166 ptr[i-n0]=s;
167 s+=itmp;
168 }
169 ptr[n-n0]=s;
170 // fill index
171 int *index=new int[ptr[n-n0]];
172 #pragma omp parallel for
173 for (int i=n0; i<n; ++i) {
174 IndexList_toArray(index_list[i], &index[ptr[i-n0]], range_min, range_max, index_offset);
175 }
176
177 Paso_Pattern* out=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, n-n0, range_max+index_offset, ptr, index);
178 if (!noError()) {
179 delete[] ptr;
180 delete[] index;
181 Paso_Pattern_free(out);
182 }
183 return out;
184 }
185
186 } // namespace finley
187

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
svn:mergeinfo /branches/lapack2681/finley/src/IndexList.cpp:2682-2741 /branches/pasowrap/finley/src/IndexList.cpp:3661-3674 /branches/py3_attempt2/finley/src/IndexList.cpp:3871-3891 /branches/restext/finley/src/IndexList.cpp:2610-2624 /branches/ripleygmg_from_3668/finley/src/IndexList.cpp:3669-3791 /branches/stage3.0/finley/src/IndexList.cpp:2569-2590 /branches/symbolic_from_3470/finley/src/IndexList.cpp:3471-3974 /release/3.0/finley/src/IndexList.cpp:2591-2601 /trunk/finley/src/IndexList.cpp:4257-4344

  ViewVC Help
Powered by ViewVC 1.1.26