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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4492 - (show annotations)
Tue Jul 2 01:44:11 2013 UTC (5 years, 9 months ago) by caltinay
File size: 9773 byte(s)
Finley changes that were held back while in release mode - moved more stuff
into finley namespace.

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 Adds points to the Points table of the mesh.
20
21 *****************************************************************************/
22
23 #include "Mesh.h"
24
25 #ifdef ESYS_MPI
26 void Finley_Mesh_MPI_minimizeDistance(void *invec_p, void *inoutvec_p,
27 int *len, MPI_Datatype *dtype)
28 {
29 const int numPoints = (*len)/2;
30 double *invec = reinterpret_cast<double*>(invec_p);
31 double *inoutvec = reinterpret_cast<double*>(inoutvec_p);
32 for (int i=0; i<numPoints; i++) {
33 if (invec[2*i] < inoutvec[2*i]) {
34 inoutvec[2*i] = invec[2*i];
35 inoutvec[2*i+1] = invec[2*i+1];
36 }
37 }
38 }
39 #endif
40
41 void Finley_Mesh_addPoints(Finley_Mesh* mesh, const dim_t numPoints,
42 const double* points_ptr, const index_t* tags_ptr)
43 {
44 Esys_MPIInfo *mpi_info = Esys_MPIInfo_getReference(mesh->MPIInfo);
45 ElementFile *oldPoints=mesh->Points;
46 ReferenceElementSet *refPoints=NULL;
47 dim_t numOldPoints;
48 if (oldPoints == NULL) {
49 refPoints=ReferenceElementSet_alloc(Point1,
50 mesh->integrationOrder, mesh->reducedIntegrationOrder);
51 numOldPoints=0;
52 } else {
53 refPoints=ReferenceElementSet_reference(oldPoints->referenceElementSet);
54 numOldPoints=mesh->Points->numElements;
55 }
56 ElementFile *newPoints=new ElementFile(refPoints, mpi_info);
57
58 // first we find the node which is the closest on this processor:
59 double *dist_p = new double[numPoints];
60 index_t *node_id_p = new index_t[numPoints];
61 index_t *point_index_p = new index_t[numPoints];
62
63 for (dim_t i=0; i<numPoints; ++i) {
64 dist_p[i]=LARGE_POSITIVE_FLOAT;
65 node_id_p[i]=-1;
66 node_id_p[i]=-1;
67 }
68
69 const double *coords = mesh->Nodes->Coordinates;
70 const dim_t numDim=mesh->Nodes->numDim;
71 if (numDim == 3) {
72 #pragma omp parallel
73 {
74 for (dim_t i=0; i<numPoints; ++i) {
75 const double X0=points_ptr[INDEX2(0,i,numDim)];
76 const double X1=points_ptr[INDEX2(1,i,numDim)];
77 const double X2=points_ptr[INDEX2(2,i,numDim)];
78 double dist_local=LARGE_POSITIVE_FLOAT;
79 index_t node_id_local=-1;
80 #pragma omp for
81 for (dim_t n=0; n<mesh->Nodes->numNodes; n++) {
82 const double D0=coords[INDEX2(0,n,numDim)] - X0;
83 const double D1=coords[INDEX2(1,n,numDim)] - X1;
84 const double D2=coords[INDEX2(2,n,numDim)] - X2;
85 const double d = D0*D0 + D1*D1 + D2*D2;
86 if (d < dist_local) {
87 dist_local = d;
88 node_id_local = n;
89 }
90 }
91 #pragma omp critical
92 {
93 if (dist_local < dist_p[i]) {
94 dist_p[i] = dist_local;
95 node_id_p[i] = node_id_local;
96 }
97 }
98 }
99 } // end parallel section
100 } else if (numDim == 2) {
101 #pragma omp parallel
102 {
103 for (dim_t i=0; i<numPoints; ++i) {
104 const double X0=points_ptr[INDEX2(0,i,numDim)];
105 const double X1=points_ptr[INDEX2(1,i,numDim)];
106 double dist_local=LARGE_POSITIVE_FLOAT;
107 index_t node_id_local=-1;
108 #pragma omp for
109 for (dim_t n=0; n<mesh->Nodes->numNodes; n++) {
110 const double D0=coords[INDEX2(0,n,numDim)] - X0;
111 const double D1=coords[INDEX2(1,n,numDim)] - X1;
112 const double d = D0*D0 + D1*D1;
113 if (d < dist_local) {
114 dist_local = d;
115 node_id_local = n;
116 }
117 }
118 #pragma omp critical
119 {
120 if (dist_local < dist_p[i]) {
121 dist_p[i] = dist_local;
122 node_id_p[i] = node_id_local;
123 }
124 }
125 }
126 } // end parallel section
127 } else { // numDim==1
128 #pragma omp parallel
129 {
130 for (dim_t i=0; i<numPoints; ++i) {
131 const double X0=points_ptr[INDEX2(0,i,numDim)];
132 double dist_local=LARGE_POSITIVE_FLOAT;
133 index_t node_id_local=-1;
134 #pragma omp for
135 for (dim_t n=0; n< mesh-> Nodes->numNodes; n++) {
136 const double D0=coords[INDEX2(0,n,numDim)] - X0;
137 const double d = D0*D0;
138 if (d < dist_local) {
139 dist_local = d;
140 node_id_local = n;
141 }
142 }
143 #pragma omp critical
144 {
145 if (dist_local < dist_p[i]) {
146 dist_p[i] = dist_local;
147 node_id_p[i] = node_id_local;
148 }
149 }
150 }
151 } // end parallel section
152 } //numDim
153
154 #ifdef ESYS_MPI
155 // now we need to reduce this across all processors
156 {
157 MPI_Op op;
158 int count = 2*numPoints;
159 double *sendbuf=new double[count];
160 double *recvbuf=new double[count];
161
162 for (dim_t i=0; i<numPoints; ++i) {
163 sendbuf[2*i ]=dist_p[i];
164 sendbuf[2*i+1]=static_cast<double>(mesh->Nodes->Id[node_id_p[i]]);
165 }
166 MPI_Op_create(Finley_Mesh_MPI_minimizeDistance, TRUE, &op);
167 MPI_Allreduce(sendbuf, recvbuf, count, MPI_DOUBLE, op, mpi_info->comm);
168 MPI_Op_free(&op);
169 // if the node id has changed we found another node which is closer
170 // elsewhere
171 for (dim_t i=0; i<numPoints; ++i) {
172 const int best_fit_Id = static_cast<const int>(recvbuf[2*i+1]+0.5);
173 if (best_fit_Id != mesh->Nodes->Id[node_id_p[i]]) {
174 node_id_p[i] = -1;
175 }
176 }
177 delete[] sendbuf;
178 delete[] recvbuf;
179 }
180 #endif
181 delete[] dist_p;
182
183 // we pick the points to be used on this processor
184 dim_t numNewPoints=0;
185 const index_t firstDOF=Paso_Distribution_getFirstComponent(mesh->Nodes->degreesOfFreedomDistribution);
186 const index_t lastDOF=Paso_Distribution_getLastComponent(mesh->Nodes->degreesOfFreedomDistribution);
187
188 for (dim_t i=0; i<numPoints; ++i) {
189 if (node_id_p[i]>-1) {
190 // this processor uses a node which is identical to point i
191 if (mesh->Nodes->globalReducedDOFIndex[node_id_p[i]] > -1) {
192 // the point is also used in the reduced mesh
193 const index_t global_id=mesh->Nodes->globalDegreesOfFreedom[node_id_p[i]];
194 if ( (firstDOF<=global_id) && (global_id<lastDOF) ) {
195 // is this point actually relevant
196 bool notAnOldPoint=true;
197 if (numOldPoints > 0) {
198 // is this point already in the Point table?
199 for (dim_t k=0; k<numOldPoints; ++k) {
200 if (global_id == oldPoints->Nodes[k]) {
201 notAnOldPoint=false;
202 break;
203 }
204 }
205 }
206 if (notAnOldPoint) {
207 // is this point unique in the new list of points?
208 bool notANewPoint=true;
209 for (dim_t k=0; k<numNewPoints; ++k) {
210 if (global_id == mesh->Nodes->globalDegreesOfFreedom[node_id_p[point_index_p[k]]]) {
211 notANewPoint=false;
212 break;
213 }
214 }
215 if (notANewPoint) {
216 point_index_p[numNewPoints]=i;
217 numNewPoints++;
218 }
219 }
220 }
221 }
222 }
223 }
224
225 // now we are ready to create the new Point table
226 newPoints->allocTable(numOldPoints+numNewPoints);
227 if (numOldPoints > 0) {
228 #pragma omp parallel for schedule(static)
229 for(dim_t n=0; n<numOldPoints; n++) {
230 newPoints->Owner[n]=oldPoints->Owner[n];
231 newPoints->Id[n] =oldPoints->Id[n];
232 newPoints->Tag[n] =oldPoints->Tag[n];
233 newPoints->Nodes[n]=oldPoints->Nodes[n];
234 newPoints->Color[n]=0;
235 }
236 }
237 #pragma omp parallel for schedule(static)
238 for(dim_t n=0; n<numNewPoints; n++) {
239 const index_t idx = point_index_p[n];
240 newPoints->Owner[numOldPoints+n]=mpi_info->rank;
241 newPoints->Id[numOldPoints+n] =0;
242 newPoints->Tag[numOldPoints+n] =tags_ptr[idx];
243 newPoints->Nodes[numOldPoints+n]=node_id_p[idx];
244 newPoints->Color[numOldPoints+n]=0;
245 }
246 newPoints->minColor=0;
247 newPoints->maxColor=0;
248
249 // all done, clean up
250 delete[] node_id_p;
251 delete[] point_index_p;
252 ReferenceElementSet_dealloc(refPoints);
253 Esys_MPIInfo_free(mpi_info);
254 if (Finley_noError()) {
255 delete oldPoints;
256 mesh->Points=newPoints;
257 } else {
258 delete newPoints;
259 }
260 }
261

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26