test_propagators_00002.cpp
1 /*---------------------------------------------------------------------------*\
2  *
3  * mimmo
4  *
5  * Copyright (C) 2015-2021 OPTIMAD engineering Srl
6  *
7  * -------------------------------------------------------------------------
8  * License
9  * This file is part of mimmo.
10  *
11  * mimmo is free software: you can redistribute it and/or modify it
12  * under the terms of the GNU Lesser General Public License v3 (LGPL)
13  * as published by the Free Software Foundation.
14  *
15  * mimmo is distributed in the hope that it will be useful, but WITHOUT
16  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
18  * License for more details.
19  *
20  * You should have received a copy of the GNU Lesser General Public License
21  * along with mimmo. If not, see <http://www.gnu.org/licenses/>.
22  *
23  \ *---------------------------------------------------------------------------*/
24 
25 #include "mimmo_propagators.hpp"
26 
27 // =================================================================================== //
42 // =================================================================================== //
43 
44 // =================================================================================== //
45 /*
46  * Create bar grid with prismatic layer + tetrahedral bulk and an intermediate ref level between them
47  * \param[out]boundary layer surface grid pidded in 3, 0 front face prismatic, 1 back bulk bar faces, 2 other.
48  * \return the volume grid
49  */
50 mimmo::MimmoSharedPointer<mimmo::MimmoObject> createTestVolumeMesh( mimmo::MimmoSharedPointer<mimmo::MimmoObject> &boundary)
51 {
52 
53  std::array<double, 3> origin({{0.0,0.0,0.0}});
54  double width(0.2), height(0.2);
55  double prismdepth(0.4), bulkdepth(0.6);
56 
57  int nw(10), nh(10), npd(20), nbd(30);
58 
59  double deltaw = width/ double(nw);
60  double deltah = height/ double(nh);
61  double deltapd = prismdepth/ double(npd);
62  double deltabd = bulkdepth/ double(nbd);
63 
64  std::vector<std::array<double,3> > verts ((nw+1)*(nh+1)*(npd+nbd+1));
65 
66  int counter = 0;
67  //prism layer verts
68  for(int k=0; k<npd; ++k){
69  for(int j=0; j<=nh; ++j){
70  for(int i=0; i<=nw; ++i){
71  verts[counter][0] = origin[0] + deltaw * i;
72  verts[counter][1] = origin[1] + deltah * j;
73  verts[counter][2] = origin[2] + deltapd * k;
74  ++counter;
75  }
76  }
77  }
78  //bulk layer verts
79  for(int k=0; k<=nbd; ++k){
80  for(int j=0; j<=nh; ++j){
81  for(int i=0; i<=nw; ++i){
82  verts[counter][0] = origin[0] + deltaw * i;
83  verts[counter][1] = origin[1] + deltah * j;
84  verts[counter][2] = origin[2] + deltapd * npd + deltabd * k;
85  ++counter;
86  }
87  }
88  }
89 
90  //create the volume mesh mimmo.
91  mimmo::MimmoSharedPointer<mimmo::MimmoObject> mesh = mimmo::MimmoSharedPointer<mimmo::MimmoObject>(new mimmo::MimmoObject(2));
92  mesh->getPatch()->reserveVertices(counter);
93 
94  //pump up the vertices
95  for(const auto & vertex : verts){
96  mesh->addVertex(vertex); //automatic id assigned to vertices.
97  }
98  //using Kuhn decomposition in tetrahedra from hexahedron, cut on diagonal for prism layer cells from hexahedron
99  mesh->getPatch()->reserveCells(2*nw*nh*npd + 6*nw*nh*nbd);
100 
101  //create connectivities for wedge elements of prism layer
102  std::vector<long> conn(6,0);
103  for(int k=0; k<npd; ++k){
104  for(int j=0; j<nh; ++j){
105  for(int i=0; i<nw; ++i){
106  //FIRST 0-2-1-4-6-7 from elemental hexa
107  conn[0] = (nw+1)*(nh+1)*k + (nw+1)*j + i;
108  conn[1] = (nw+1)*(nh+1)*k + (nw+1)*(j+1) + i+1;
109  conn[2] = (nw+1)*(nh+1)*k + (nw+1)*j + i+1;
110  conn[3] = (nw+1)*(nh+1)*(k+1) + (nw+1)*j + i;
111  conn[4] = (nw+1)*(nh+1)*(k+1) + (nw+1)*(j+1) + i+1;
112  conn[5] = (nw+1)*(nh+1)*(k+1) + (nw+1)*j + i+1;
113 
114  mesh->addConnectedCell(conn, bitpit::ElementType::WEDGE);
115  //SECOND 0-3-2-4-7-6 from elemental hexa
116  conn[0] = (nw+1)*(nh+1)*k + (nw+1)*j + i;
117  conn[1] = (nw+1)*(nh+1)*k + (nw+1)*(j+1) + i;
118  conn[2] = (nw+1)*(nh+1)*k + (nw+1)*(j+1) + i+1;
119  conn[3] = (nw+1)*(nh+1)*(k+1) + (nw+1)*j + i;
120  conn[4] = (nw+1)*(nh+1)*(k+1) + (nw+1)*(j+1) + i;
121  conn[5] = (nw+1)*(nh+1)*(k+1) + (nw+1)*(j+1) + i+1;
122 
123  mesh->addConnectedCell(conn, bitpit::ElementType::WEDGE);
124  }
125  }
126  }
127 
128  //create connectivities for tetra elements: USING KUHN
129  //BEWARE, the coarser the tetrahedral kuhn adapted part the worse the non-orthogonality of cells.
130  // This impacts on LAPLACIAL SOLUTION. //TODO Use a better tetrahedral element generation.
131 
132  conn.clear();
133  conn.resize(4,0);
134  for(int k=npd; k<(npd+nbd); ++k){
135  for(int j=0; j<nh; ++j){
136  for(int i=0; i<nw; ++i){
137  //FIRST 0-1-2-6 from elemenhal hexa
138  conn[0] = (nw+1)*(nh+1)*k + (nw+1)*j + i;
139  conn[1] = (nw+1)*(nh+1)*k + (nw+1)*j + i+1;
140  conn[2] = (nw+1)*(nh+1)*k + (nw+1)*(j+1) + i+1;
141  conn[3] = (nw+1)*(nh+1)*(k+1) + (nw+1)*(j+1) + i+1;
142 
143  mesh->addConnectedCell(conn, bitpit::ElementType::TETRA);
144 
145  //SECOND 0-5-1-6 from elemenhal hexa
146  conn[0] = (nw+1)*(nh+1)*k + (nw+1)*j + i;
147  conn[1] = (nw+1)*(nh+1)*(k+1) + (nw+1)*j + i+1;
148  conn[2] = (nw+1)*(nh+1)*k + (nw+1)*j + i+1;
149  conn[3] = (nw+1)*(nh+1)*(k+1) + (nw+1)*(j+1) + i+1;
150 
151  mesh->addConnectedCell(conn, bitpit::ElementType::TETRA);
152 
153  //THIRD 0-2-3-6 from elemenhal hexa
154  conn[0] = (nw+1)*(nh+1)*k + (nw+1)*j + i;
155  conn[1] = (nw+1)*(nh+1)*k + (nw+1)*(j+1) + i+1;
156  conn[2] = (nw+1)*(nh+1)*k + (nw+1)*(j+1) + i;
157  conn[3] = (nw+1)*(nh+1)*(k+1) + (nw+1)*(j+1) + i+1;
158 
159  mesh->addConnectedCell(conn, bitpit::ElementType::TETRA);
160 
161  //FIRST 0-4-5-6 from elemenhal hexa
162  conn[0] = (nw+1)*(nh+1)*k + (nw+1)*j + i;
163  conn[1] = (nw+1)*(nh+1)*(k+1) + (nw+1)*j + i;
164  conn[2] = (nw+1)*(nh+1)*(k+1) + (nw+1)*j + i+1;
165  conn[3] = (nw+1)*(nh+1)*(k+1) + (nw+1)*(j+1) + i+1;
166 
167  mesh->addConnectedCell(conn, bitpit::ElementType::TETRA);
168 
169  //SECOND 0-7-4-6 from elemenhal hexa
170  conn[0] = (nw+1)*(nh+1)*k + (nw+1)*j + i;
171  conn[1] = (nw+1)*(nh+1)*(k+1) + (nw+1)*(j+1) + i;
172  conn[2] = (nw+1)*(nh+1)*(k+1) + (nw+1)*j + i;
173  conn[3] = (nw+1)*(nh+1)*(k+1) + (nw+1)*(j+1) + i+1;
174 
175  mesh->addConnectedCell(conn, bitpit::ElementType::TETRA);
176 
177  //THIRD 0-3-7-6 from elemenhal hexa
178  conn[0] = (nw+1)*(nh+1)*k + (nw+1)*j + i;
179  conn[1] = (nw+1)*(nh+1)*k + (nw+1)*(j+1) + i;
180  conn[2] = (nw+1)*(nh+1)*(k+1) + (nw+1)*(j+1) + i;
181  conn[3] = (nw+1)*(nh+1)*(k+1) + (nw+1)*(j+1) + i+1;
182 
183  mesh->addConnectedCell(conn, bitpit::ElementType::TETRA);
184 
185  }
186  }
187  }
188 
189  mesh->updateAdjacencies();
190  mesh->updateInterfaces();
191  mesh->update();
192 
193  livector1D borderverts = mesh->extractBoundaryVertexID();
194 
195  //get the subset of ids nearest at z=0 and z= prismdepth+bulkdepth;
196  double zin(0.0), zout(prismdepth+bulkdepth), zwork;
197  std::vector<long> boundary1verts, boundary2verts;
198  for(long idv : borderverts){
199  zwork = mesh->getVertexCoords(idv)[2];
200  if(std::abs(zwork - zin) < 0.5*deltapd){
201  boundary1verts.push_back(idv);
202  }
203  if(std::abs(zwork - zout) < 0.5*deltabd){
204  boundary2verts.push_back(idv);
205  }
206  }
207  //serial test
208  livector1D boundary1faces = mesh->getInterfaceFromVertexList(boundary1verts, true, true);
209  livector1D boundary2faces = mesh->getInterfaceFromVertexList(boundary2verts, true, true);
210 
211  //put together
212  mimmo::MimmoPiercedVector<long> pidfaces;
213  pidfaces.reserve(boundary1faces.size()+boundary2faces.size());
214  for(long id : boundary1faces){
215  pidfaces.insert(id,0);
216  }
217  for(long id : boundary2faces){
218  pidfaces.insert(id,1);
219  }
220 
221  std::vector<long> boundaryverts(boundary1verts.begin(), boundary1verts.end());
222  boundaryverts.insert(boundaryverts.end(), boundary2verts.begin(), boundary2verts.end());
223 
224  //create the boundary meshes
225  auto orinterfaces = mesh->getInterfaces();
226  auto orcells = mesh->getCells();
227 
228  boundary = mimmo::MimmoSharedPointer<mimmo::MimmoObject>(new mimmo::MimmoObject(1));
229  boundary->getPatch()->reserveVertices(boundaryverts.size());
230  boundary->getPatch()->reserveCells(pidfaces.size());
231 
232  //push in verts
233  for(long id : boundaryverts){
234  boundary->addVertex(mesh->getVertexCoords(id), id);
235  }
236 
237  //push in interfaces ad 2D cells
238  for(auto it=pidfaces.begin(); it!=pidfaces.end(); ++it){
239  bitpit::Interface & ii = orinterfaces.at(it.getId());
240  long * conn = ii.getConnect();
241  std::size_t connsize = ii.getConnectSize();
242 
243  bitpit::ElementType et = orcells.at(ii.getOwner()).getFaceType(ii.getOwnerFace());
244 
245  boundary->addConnectedCell(std::vector<long>(&conn[0], &conn[connsize]), et, *it,it.getId());
246  }
247  boundary->updateAdjacencies();
248  boundary->update();
249 
250  return mesh;
251 }
252 
253 
254 // =================================================================================== //
255 /*
256  Testing scalar field propagation ond mixed element and messy mesh.
257 */
258 
259 int test1() {
260 
261  mimmo::MimmoSharedPointer<mimmo::MimmoObject> boundary;
262  mimmo::MimmoSharedPointer<mimmo::MimmoObject> mesh = createTestVolumeMesh(boundary);
263 
264  bool check = false;
265  long targetNode = (11*11)*25 + (11)*5 + 5;
266 
267 //TESTING THE SCALAR PROPAGATOR //////
268  // and the scalar field of Dirichlet values on its nodes.
269  mimmo::MimmoPiercedVector<double> bc_surf_field;
270  bc_surf_field.setGeometry(boundary);
271  bc_surf_field.setDataLocation(mimmo::MPVLocation::POINT);
272  bc_surf_field.reserve(boundary->getNVertices());
273  for(long val : boundary->getVertices().getIds()){
274  bc_surf_field.insert(val, 0.0);
275  }
276 
277  //extract vertex of boundary cells pidded 0 and assign to them a condition 10.0
278  std::vector<long> listPP = boundary->getVertexFromCellList(boundary->extractPIDCells(0));
279  for(long id : listPP){
280  bc_surf_field.at(id) = 10.0;
281  }
282 
283  // Now create a PropagateScalarField and solve the laplacian.
284  mimmo::PropagateScalarField * prop = new mimmo::PropagateScalarField();
285  prop->setName("test00002_PropagateScalarField");
286  prop->setGeometry(mesh);
287  prop->addDirichletBoundaryPatch(boundary);
288  prop->addDirichletConditions(&bc_surf_field);
289  prop->setDamping(false);
290  prop->setPlotInExecution(true);
291  prop->setSolverMultiStep(4);
292  prop->exec();
293 
294  auto values = prop->getPropagatedField();
295  std::cout<< values->at(targetNode) <<std::endl;
296  check = check || (std::abs(values->at(targetNode)-3.75446) > 1.0E-5);
297 // check = false;
298 
299  delete prop;
300 
301  return check;
302 }
303 
304 // =================================================================================== //
305 
306 int main( int argc, char *argv[] ) {
307 
308  BITPIT_UNUSED(argc);
309  BITPIT_UNUSED(argv);
310 
311 #if MIMMO_ENABLE_MPI
312  MPI_Init(&argc, &argv);
313 #endif
314 
315  int val = 1;
316  try{
317  val = test1() ;
318  }
319  catch(std::exception & e){
320  std::cout<<"test_propagators_00002 exited with an error of type : "<<e.what()<<std::endl;
321  return 1;
322  }
323 #if MIMMO_ENABLE_MPI
324  MPI_Finalize();
325 #endif
326 
327  return val;
328 }
mimmo::MimmoObject
MimmoObject is the basic geometry container for mimmo library.
Definition: MimmoObject.hpp:143
livector1D
std::vector< long > livector1D
Definition: mimmoTypeDef.hpp:42
mimmo::MimmoPiercedVector::setGeometry
void setGeometry(MimmoSharedPointer< MimmoObject > geo)
Definition: MimmoPiercedVector.tpp:314
mimmo::MPVLocation::POINT
@ POINT
mimmo::MimmoPiercedVector::setDataLocation
void setDataLocation(MPVLocation loc)
Definition: MimmoPiercedVector.tpp:324
Generated on Tue Mar 16 2021 16:26:41 for mimmo by   doxygen 1.8.17