Line data Source code
1 : /*
2 : * Copyright (c) 2015: G-CSC, Goethe University Frankfurt
3 : * Author: Sebastian Reiter
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #include "field_util.h"
34 : #include "lib_grid/algorithms/geom_obj_util/vertex_util.h"
35 :
36 : using namespace std;
37 :
38 : namespace ug{
39 :
40 : ////////////////////////////////////////////////////////////////////////////////
41 0 : void CreateGridFromField(Grid& grid,
42 : const Field<number>& field,
43 : const vector2& cellSize,
44 : const vector2& offset,
45 : number noDataValue,
46 : Grid::VertexAttachmentAccessor<APosition> aaPos)
47 : {
48 0 : Field<Vertex*> vrtField(field.width(), field.height());
49 : vrtField.fill_all(0);
50 :
51 : // iterate over cells and create triangles or quadrilaterals between
52 : // adjacent vertices
53 0 : const int fieldWidth = (int)field.width();
54 0 : const int fieldHeight = (int)field.height();
55 :
56 0 : for(int irow = 0; irow + 1 < fieldHeight; ++irow){
57 0 : for(int icol = 0; icol + 1 < fieldWidth; ++icol){
58 : // corner indices in ccw order
59 : // index order:
60 : // 0--3
61 : // | |
62 : // 1--2
63 : const int ci[4][2] = {{icol, irow},
64 : {icol, irow + 1},
65 : {icol + 1, irow + 1},
66 0 : {icol + 1, irow}};
67 :
68 : // cell corner values in ccw order
69 : number val[4];
70 : size_t numVals = 0;
71 : int firstInvalid = -1;
72 0 : for(size_t i = 0; i < 4; ++i){
73 0 : val[i] = field.at(ci[i][0], ci[i][1]);
74 0 : if(val[i] != noDataValue)
75 0 : ++numVals;
76 0 : else if(firstInvalid == -1)
77 0 : firstInvalid = (int)i;
78 : }
79 :
80 0 : if(numVals < 3)
81 0 : continue;
82 :
83 : // create vertices as necessary and store corner vertices in vrts[]
84 : Vertex* vrts[4];
85 0 : for(size_t i = 0; i < 4; ++i){
86 0 : const int ic = ci[i][0];
87 0 : const int ir = ci[i][1];
88 0 : Vertex* vrt = vrtField.at(ic, ir);
89 0 : if((val[i] != noDataValue) && (!vrt)){
90 0 : vrt = *grid.create<RegularVertex>();
91 : aaPos[vrt] =
92 0 : vector3(offset.x() + cellSize.x() * (number)ic,
93 0 : offset.y() + cellSize.y() * (number)ir,
94 : field.at(ic, ir));
95 0 : vrtField.at(ic, ir) = vrt;
96 : }
97 0 : vrts[i] = vrt;
98 : }
99 :
100 0 : if(numVals == 4){
101 : // create a quad
102 0 : grid.create<Quadrilateral>(
103 0 : QuadrilateralDescriptor(vrts[3], vrts[2], vrts[1], vrts[0]));
104 : }
105 0 : else if(numVals == 3){
106 : // create a tri
107 0 : grid.create<Triangle>(
108 0 : TriangleDescriptor(vrts[(firstInvalid + 3) % 4],
109 0 : vrts[(firstInvalid + 2) % 4],
110 0 : vrts[(firstInvalid + 1) % 4]));
111 : }
112 : }
113 : }
114 0 : }
115 :
116 : ////////////////////////////////////////////////////////////////////////////////
117 0 : void CreateGridFromFieldBoundary(Grid& grid,
118 : const Field<number>& field,
119 : const vector2& cellSize,
120 : const vector2& offset,
121 : number noDataValue,
122 : Grid::VertexAttachmentAccessor<APosition> aaPos)
123 : {
124 0 : Field<Vertex*> vrtField(field.width(), field.height());
125 : vrtField.fill_all(0);
126 :
127 : // iterate over cells and create triangles or quadrilaterals between
128 : // adjacent vertices
129 0 : const int fieldWidth = (int)field.width();
130 0 : const int fieldHeight = (int)field.height();
131 :
132 0 : for(int irow = 0; irow + 1 < fieldHeight; ++irow){
133 0 : for(int icol = 0; icol + 1 < fieldWidth; ++icol){
134 : // corner indices in ccw order
135 : // index order:
136 : // 0--3
137 : // | |
138 : // 1--2
139 : const int ci[4][2] = {{icol, irow},
140 : {icol, irow + 1},
141 : {icol + 1, irow + 1},
142 0 : {icol + 1, irow}};
143 :
144 : // cell corner values in ccw order
145 : number val[4];
146 : size_t numVals = 0;
147 : int firstInvalid = -1;
148 0 : for(size_t i = 0; i < 4; ++i){
149 0 : val[i] = field.at(ci[i][0], ci[i][1]);
150 0 : if(val[i] != noDataValue)
151 0 : ++numVals;
152 0 : else if(firstInvalid == -1)
153 0 : firstInvalid = (int)i;
154 : }
155 :
156 0 : if(numVals < 3)
157 0 : continue;
158 :
159 : // create vertices as necessary and store corner vertices in vrts[]
160 : Vertex* vrts[4];
161 :
162 0 : if(numVals == 3){
163 : // create a diagonal which doesn't include the firstInvalid corner
164 0 : int ind[2] = {(firstInvalid + 1) % 4, (firstInvalid + 3) % 4};
165 :
166 0 : for(int i = 0; i < 2; ++i){
167 0 : int ic = ci[ind[i]][0];
168 0 : int ir = ci[ind[i]][1];
169 0 : Vertex* vrt = vrtField.at(ic, ir);
170 0 : if((val[ind[i]] != noDataValue) && (!vrt)){
171 0 : vrt = *grid.create<RegularVertex>();
172 : aaPos[vrt] =
173 0 : vector3(offset.x() + cellSize.x() * (number)ic,
174 0 : offset.y() + cellSize.y() * (number)ir,
175 : val[ind[i]]);
176 0 : vrtField.at(ic, ir) = vrt;
177 : }
178 0 : vrts[ind[i]] = vrt;
179 : }
180 :
181 0 : grid.create<RegularEdge>(EdgeDescriptor(vrts[ind[0]], vrts[ind[1]]));
182 : }
183 :
184 : // now check for each inner side whether it's a neighbor to an empty cell
185 0 : const int colAdd[4] = {-1, 0, 1, 0};
186 0 : const int rowAdd[4] = {0, 1, 0, -1};
187 0 : for(int iside = 0; iside < 4; ++iside){
188 0 : const int ind[2] = {iside, (iside + 1) % 4};
189 0 : if(val[ind[0]] != noDataValue && val[ind[1]] != noDataValue){
190 : // the edge is either inner or boundary
191 : bool gotOne = false;
192 0 : for(int i = 0; i < 2; ++i){
193 0 : const int ic = ci[ind[i]][0] + colAdd[iside];
194 0 : const int ir = ci[ind[i]][1] + rowAdd[iside];
195 0 : if(ic >= 0 && ir >= 0 && ic < fieldWidth && ir < fieldHeight){
196 0 : if(field.at(ic, ir) != noDataValue){
197 : gotOne = true;
198 : break;
199 : }
200 : }
201 : }
202 :
203 0 : if(!gotOne){
204 : // the current side has to be represented by an edge
205 0 : for(int i = 0; i < 2; ++i){
206 0 : const int ic = ci[ind[i]][0];
207 0 : const int ir = ci[ind[i]][1];
208 0 : Vertex* vrt = vrtField.at(ic, ir);
209 0 : if((val[ind[i]] != noDataValue) && (!vrt)){
210 0 : vrt = *grid.create<RegularVertex>();
211 : aaPos[vrt] =
212 0 : vector3(offset.x() + cellSize.x() * (number)ic,
213 0 : offset.y() + cellSize.y() * (number)ir,
214 : val[ind[i]]);
215 0 : vrtField.at(ic, ir) = vrt;
216 : }
217 0 : vrts[ind[i]] = vrt;
218 : }
219 0 : grid.create<RegularEdge>(EdgeDescriptor(vrts[ind[0]], vrts[ind[1]]));
220 : }
221 : }
222 : }
223 : }
224 : }
225 0 : }
226 :
227 : }// end of namespace
|
