/[escript]/branches/trilinos_from_5897/dudley/src/Util.cpp
ViewVC logotype

Contents of /branches/trilinos_from_5897/dudley/src/Util.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6079 - (show annotations)
Mon Mar 21 12:22:38 2016 UTC (2 years, 10 months ago) by caltinay
File size: 7929 byte(s)
Big commit - making dudley much more like finley to make it more
managable. Fixed quite a few issues that had been fixed in finley.
Disposed of all ReducedNode/ReducedDOF entities that dudley never supported.
Compiles and passes tests.

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2016 by The 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 2012-2013 by School of Earth Sciences
13 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 *
15 *****************************************************************************/
16
17 #include "Util.h"
18
19 #include <escript/index.h>
20
21 #include <algorithm> // std::sort
22
23 namespace dudley {
24 namespace util {
25
26 /// comparison function for sortValueAndIndex
27 bool ValueAndIndexCompare(const std::pair<int,int> &i, const std::pair<int, int> &j)
28 {
29 // to ensure we have a strict ordering as required by std
30 if (i.first == j.first)
31 return i.second < j.second;
32 return i.first < j.first;
33 }
34
35 void sortValueAndIndex(ValueAndIndexList& array)
36 {
37 std::sort(array.begin(), array.end(), ValueAndIndexCompare);
38 }
39
40 void gather(int len, const index_t* index, int numData, const double* in,
41 double* out)
42 {
43 for (int s = 0; s < len; s++) {
44 for (int i = 0; i < numData; i++) {
45 out[INDEX2(i, s, numData)] = in[INDEX2(i, index[s], numData)];
46 }
47 }
48 }
49
50 void addScatter(int len, const index_t* index, int numData,
51 const double* in, double* out, index_t upperBound)
52 {
53 for (int s = 0; s < len; s++) {
54 for (int i = 0; i < numData; i++) {
55 if (index[s] < upperBound) {
56 out[INDEX2(i, index[s], numData)] += in[INDEX2(i, s, numData)];
57 }
58 }
59 }
60 }
61
62 void smallMatMult(int A1, int A2, double* A, int B2, const double* B,
63 const double* C)
64 {
65 for (int i = 0; i < A1; i++) {
66 for (int j = 0; j < A2; j++) {
67 double sum = 0.;
68 for (int s = 0; s < B2; s++)
69 sum += B[INDEX2(i,s,A1)] * C[INDEX2(s,j,B2)];
70 A[INDEX2(i,j,A1)] = sum;
71 }
72 }
73 }
74
75 void smallMatSetMult1(int len, int A1, int A2, double* A, int B2,
76 const double* B, const double* C)
77 {
78 for (int q = 0; q < len; q++) {
79 for (int i = 0; i < A1; i++) {
80 for (int j = 0; j < A2; j++) {
81 double sum = 0.;
82 for (int s = 0; s < B2; s++)
83 sum += B[INDEX3(i,s,q,A1,B2)] * C[INDEX2(s,j,B2)];
84 A[INDEX3(i,j,q,A1,A2)] = sum;
85 }
86 }
87 }
88 }
89
90 void normalVector(int len, int dim, int dim1, const double* A, double* Normal)
91 {
92 int q;
93
94 switch (dim) {
95 case 1:
96 for (q = 0; q < len; q++)
97 Normal[q] = 1.;
98 break;
99 case 2:
100 for (q = 0; q < len; q++) {
101 const double A11 = A[INDEX3(0,0,q,2,dim1)];
102 const double A21 = A[INDEX3(1,0,q,2,dim1)];
103 const double length = sqrt(A11*A11+A21*A21);
104 if (length <= 0) {
105 throw DudleyException("normalVector: area equals zero.");
106 } else {
107 const double invlength = 1./length;
108 Normal[INDEX2(0,q,2)] = A21*invlength;
109 Normal[INDEX2(1,q,2)] = -A11*invlength;
110 }
111 }
112 break;
113 case 3:
114 for (q = 0; q < len; q++) {
115 const double A11 = A[INDEX3(0,0,q,3,dim1)];
116 const double A21 = A[INDEX3(1,0,q,3,dim1)];
117 const double A31 = A[INDEX3(2,0,q,3,dim1)];
118 const double A12 = A[INDEX3(0,1,q,3,dim1)];
119 const double A22 = A[INDEX3(1,1,q,3,dim1)];
120 const double A32 = A[INDEX3(2,1,q,3,dim1)];
121 const double CO_A13 = A21*A32-A31*A22;
122 const double CO_A23 = A31*A12-A11*A32;
123 const double CO_A33 = A11*A22-A21*A12;
124 const double length = sqrt(CO_A13*CO_A13 + CO_A23*CO_A23
125 + CO_A33*CO_A33);
126 if (length <= 0) {
127 throw DudleyException("normalVector: area equals zero.");
128 } else {
129 const double invlength = 1./length;
130 Normal[INDEX2(0,q,3)] = CO_A13*invlength;
131 Normal[INDEX2(1,q,3)] = CO_A23*invlength;
132 Normal[INDEX2(2,q,3)] = CO_A33*invlength;
133 }
134 }
135 break;
136 }
137 }
138
139 IndexPair getMinMaxInt(int dim, dim_t N, const index_t* values)
140 {
141 index_t vmin = escript::DataTypes::index_t_max();
142 index_t vmax = escript::DataTypes::index_t_min();
143 if (values && dim*N > 0) {
144 vmin = vmax = values[0];
145 #pragma omp parallel
146 {
147 index_t vmin_local = vmin;
148 index_t vmax_local = vmax;
149 #pragma omp for
150 for (index_t j = 0; j < N; j++) {
151 for (int i = 0; i < dim; i++) {
152 vmin_local = std::min(vmin_local, values[INDEX2(i,j,dim)]);
153 vmax_local = std::max(vmax_local, values[INDEX2(i,j,dim)]);
154 }
155 }
156 #pragma omp critical
157 {
158 vmin = std::min(vmin_local, vmin);
159 vmax = std::max(vmax_local, vmax);
160 }
161 }
162 }
163 return IndexPair(vmin,vmax);
164 }
165
166 IndexPair getFlaggedMinMaxInt(dim_t N, const index_t* values, index_t ignore)
167 {
168 index_t vmin = escript::DataTypes::index_t_max();
169 index_t vmax = escript::DataTypes::index_t_min();
170 if (values && N > 0) {
171 vmin = vmax = values[0];
172 #pragma omp parallel
173 {
174 index_t vmin_local = vmin;
175 index_t vmax_local = vmax;
176 #pragma omp for
177 for (index_t i = 0; i < N; i++) {
178 if (values[i] != ignore) {
179 vmin_local = std::min(vmin_local, values[i]);
180 vmax_local = std::max(vmax_local, values[i]);
181 }
182 }
183 #pragma omp critical
184 {
185 vmin = std::min(vmin_local, vmin);
186 vmax = std::max(vmax_local, vmax);
187 }
188 }
189 }
190 return IndexPair(vmin,vmax);
191 }
192
193 std::vector<index_t> packMask(const std::vector<short>& mask)
194 {
195 std::vector<index_t> index;
196 for (index_t k = 0; k < mask.size(); k++) {
197 if (mask[k] >= 0) {
198 index.push_back(k);
199 }
200 }
201 return index;
202 }
203
204 void setValuesInUse(const int* values, dim_t numValues,
205 std::vector<int>& valuesInUse, escript::JMPI mpiinfo)
206 {
207 const int MAX_VALUE = std::numeric_limits<int>::max();
208 int lastFoundValue = std::numeric_limits<int>::min();
209 bool allFound = false;
210
211 valuesInUse.clear();
212
213 while (!allFound) {
214 // find smallest value bigger than lastFoundValue
215 int minFoundValue = MAX_VALUE;
216 #pragma omp parallel
217 {
218 int local_minFoundValue = minFoundValue;
219 #pragma omp for
220 for (index_t i = 0; i < numValues; i++) {
221 const int val = values[i];
222 if (val > lastFoundValue && val < local_minFoundValue)
223 local_minFoundValue = val;
224 }
225 #pragma omp critical
226 {
227 if (local_minFoundValue < minFoundValue)
228 minFoundValue = local_minFoundValue;
229 }
230 }
231 #ifdef ESYS_MPI
232 int local_minFoundValue = minFoundValue;
233 MPI_Allreduce(&local_minFoundValue, &minFoundValue, 1, MPI_INT,
234 MPI_MIN, mpiinfo->comm);
235 #endif
236
237 // if we found a new value we need to add this to valuesInUse
238 if (minFoundValue < MAX_VALUE) {
239 valuesInUse.push_back(minFoundValue);
240 lastFoundValue = minFoundValue;
241 } else {
242 allFound = true;
243 }
244 }
245 }
246
247 } // namespace util
248 } // namespace dudley
249

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
svn:mergeinfo /branches/4.0fordebian/dudley/src/Util.cpp:5567-5588 /branches/complex/dudley/src/Util.cpp:5866-5937 /branches/lapack2681/finley/src/Util.cpp:2682-2741 /branches/pasowrap/dudley/src/Util.cpp:3661-3674 /branches/py3_attempt2/dudley/src/Util.cpp:3871-3891 /branches/restext/finley/src/Util.cpp:2610-2624 /branches/ripleygmg_from_3668/dudley/src/Util.cpp:3669-3791 /branches/stage3.0/finley/src/Util.cpp:2569-2590 /branches/symbolic_from_3470/dudley/src/Util.cpp:3471-3974 /branches/symbolic_from_3470/ripley/test/python/dudley/src/Util.cpp:3517-3974 /release/3.0/finley/src/Util.cpp:2591-2601 /release/4.0/dudley/src/Util.cpp:5380-5406 /trunk/dudley/src/Util.cpp:4257-4344,5898-6007 /trunk/ripley/test/python/dudley/src/Util.cpp:3480-3515

  ViewVC Help
Powered by ViewVC 1.1.26