1 /*****************************************************************************
2 * *
3 * Elmer, A Finite Element Software for Multiphysical Problems *
4 * *
5 * Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland *
6 * *
7 * This program is free software; you can redistribute it and/or *
8 * modify it under the terms of the GNU General Public License *
9 * as published by the Free Software Foundation; either version 2 *
10 * of the License, or (at your option) any later version. *
11 * *
12 * This program is distributed in the hope that it will be useful, *
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15 * GNU General Public License for more details. *
16 * *
17 * You should have received a copy of the GNU General Public License *
18 * along with this program (in file fem/GPL-2); if not, write to the *
19 * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, *
20 * Boston, MA 02110-1301, USA. *
21 * *
22 *****************************************************************************/
23
24 /*****************************************************************************
25 * *
26 * ELMER/Mesh3D nglib_api *
27 * *
28 *****************************************************************************
29 * *
30 * Authors: Mikko Lyly, Juha Ruokolainen and Peter Raback *
31 * Email: Juha.Ruokolainen@csc.fi *
32 * Web: http://www.csc.fi/elmer *
33 * Address: CSC - IT Center for Science Ltd. *
34 * Keilaranta 14 *
35 * 02101 Espoo, Finland *
36 * *
37 * Original Date: 15 Mar 2008 *
38 * *
39 *****************************************************************************/
40
41 #include <iostream>
42 #include "nglib_api.h"
43
44 using namespace std;
45
NglibAPI()46 NglibAPI::NglibAPI()
47 {
48 }
49
50
~NglibAPI()51 NglibAPI::~NglibAPI()
52 {
53 }
54
setDim(int ngDim)55 void NglibAPI::setDim(int ngDim)
56 {
57 this->ngDim = ngDim;
58 }
59
getDim() const60 int NglibAPI::getDim() const
61 {
62 return this->ngDim;
63 }
64
setNgmesh(nglib::Ng_Mesh * ngmesh)65 void NglibAPI::setNgmesh(nglib::Ng_Mesh *ngmesh)
66 {
67 this->ngmesh = ngmesh;
68 }
69
setNggeom2D(nglib::Ng_Geometry_2D * geom2d)70 void NglibAPI::setNggeom2D(nglib::Ng_Geometry_2D *geom2d)
71 {
72 this->geom2d = geom2d;
73 }
74
75 // Populate elmer's mesh structure:
76 //-----------------------------------------------------------------------------
createElmerMeshStructure()77 mesh_t* NglibAPI::createElmerMeshStructure()
78 {
79 mesh_t *mesh = new mesh_t;
80
81 mesh->setNodes(0);
82 mesh->setPoints(0);
83 mesh->setEdges(0);
84 mesh->setSurfaces(0);
85 mesh->setElements(0);
86
87 bool twod = (ngDim == 2) ? true : false;
88
89 if(twod) {
90 create2D(mesh);
91 } else {
92 create3D(mesh);
93 }
94
95 return mesh;
96 }
97
create2D(mesh_t * mesh)98 void NglibAPI::create2D(mesh_t *mesh)
99 {
100 Meshutils meshutils;
101
102 // Node points:
103 //--------------
104 mesh->setNodes(nglib::Ng_GetNP_2D(ngmesh));
105 mesh->newNodeArray(mesh->getNodes());
106
107 for(int i = 0; i < mesh->getNodes(); i++) {
108 node_t *node = mesh->getNode(i);
109
110 double *x = node->getXvec();
111 x[0] = 0; x[1] = 0; x[2] = 0;
112 nglib::Ng_GetPoint_2D(ngmesh, i + 1, x);
113
114 node->setIndex(-1); // default
115 }
116
117 // Boundary elements:
118 //--------------------
119 mesh->setEdges(nglib::Ng_GetNSeg_2D(ngmesh));
120 mesh->newEdgeArray(mesh->getEdges());
121
122 for(int i = 0; i < mesh->getEdges(); i++) {
123 edge_t *edge = mesh->getEdge(i);
124
125 edge->setNature(PDE_BOUNDARY);
126 edge->setCode(202);
127 edge->setNodes(2);
128 edge->newNodeIndexes(2);
129 edge->setPoints(2);
130 edge->newPointIndexes(2);
131
132 edge->setPointIndex(0, -1);
133 edge->setPointIndex(1, -1);
134
135 int matIdx;
136 nglib::Ng_GetSegment_2D(ngmesh, i + 1, edge->getNodeIndexes(), &matIdx);
137
138 int bcIdx;
139 nglib::EG_GetSegmentBCProperty(ngmesh, geom2d, matIdx-1, &bcIdx);
140
141 edge->setIndex(bcIdx);
142
143 edge->setNodeIndex(0, edge->getNodeIndex(0) - 1);
144 edge->setNodeIndex(1, edge->getNodeIndex(1) - 1);
145
146 // swap orientation:
147 //------------------
148 int tmp = edge->getNodeIndex(0);
149 edge->setNodeIndex(0, edge->getNodeIndex(1));
150 edge->setNodeIndex(1, tmp);
151 }
152
153 // Elements:
154 //-----------
155 mesh->setSurfaces(nglib::Ng_GetNE_2D(ngmesh));
156 mesh->newSurfaceArray(mesh->getSurfaces());
157
158 double n[3];
159 n[0] = 0; n[1] = 0; n[2] = -1;
160
161 for(int i = 0; i < mesh->getSurfaces(); i++) {
162 surface_t *surface = mesh->getSurface(i);
163
164 surface->setNature(PDE_BULK);
165 surface->setCode(303);
166 surface->setNodes(3);
167 surface->newNodeIndexes(3);
168
169 int matIdx;
170
171 nglib::Ng_GetElement_2D(ngmesh, i+1, surface->getNodeIndexes(), &matIdx);
172
173 surface->setIndex(matIdx);
174
175 surface->setNodeIndex(0, surface->getNodeIndex(0) - 1);
176 surface->setNodeIndex(1, surface->getNodeIndex(1) - 1);
177 surface->setNodeIndex(2, surface->getNodeIndex(2) - 1);
178
179 surface->setNormalVec(n);
180 }
181
182 // Find parents for edge elements:
183 //---------------------------------
184 meshutils.findEdgeElementParents(mesh);
185
186 mesh->setDim(ngDim);
187 mesh->setCdim(ngDim);
188 }
189
create3D(mesh_t * mesh)190 void NglibAPI::create3D(mesh_t *mesh)
191 {
192 Meshutils meshutils;
193
194 // Node points:
195 //--------------
196 mesh->setNodes(nglib::Ng_GetNP(ngmesh));
197 mesh->newNodeArray(mesh->getNodes());
198
199 for(int i = 0; i < mesh->getNodes(); i++) {
200 node_t *node = mesh->getNode(i);
201 nglib::Ng_GetPoint(ngmesh, i+1, node->getXvec());
202 node->setIndex(-1); // default
203 }
204
205 // Boundary elements:
206 //--------------------
207 mesh->setSurfaces(nglib::Ng_GetNSE(ngmesh));
208 mesh->newSurfaceArray(mesh->getSurfaces());
209
210 for(int i = 0; i < mesh->getSurfaces(); i++) {
211 surface_t *surface = mesh->getSurface(i);
212
213 surface->setNature(PDE_BOUNDARY);
214 surface->setCode(303);
215 surface->setNodes(3);
216 surface->newNodeIndexes(3);
217 surface->setEdges(3);
218 surface->newEdgeIndexes(3);
219
220 int face = nglib::EG_GetSurfaceElementBCProperty(ngmesh, i+1);
221
222 surface->setIndex(face);
223
224 surface->setEdgeIndex(0, -1);
225 surface->setEdgeIndex(1, -1);
226 surface->setEdgeIndex(2, -1);
227
228 nglib::Ng_GetSurfaceElement(ngmesh, i+1, surface->getNodeIndexes());
229
230 surface->setNodeIndex(0, surface->getNodeIndex(0) - 1);
231 surface->setNodeIndex(1, surface->getNodeIndex(1) - 1);
232 surface->setNodeIndex(2, surface->getNodeIndex(2) - 1);
233
234 // swap orientation:
235 //------------------
236 int tmp = surface->getNodeIndex(1);
237 surface->setNodeIndex(1, surface->getNodeIndex(2));
238 surface->setNodeIndex(2, tmp);
239 }
240
241 // Elements:
242 //-----------
243 mesh->setElements(nglib::Ng_GetNE(ngmesh));
244 mesh->newElementArray(mesh->getElements());
245
246 for(int i = 0; i < mesh->getElements(); i++) {
247 element_t *element = mesh->getElement(i);
248
249 element->setNature(PDE_BULK);
250 element->setCode(504);
251 element->setNodes(4);
252 element->newNodeIndexes(4);
253
254 nglib::Ng_GetVolumeElement(ngmesh, i+1, element->getNodeIndexes());
255
256 element->setNodeIndex(0, element->getNodeIndex(0) - 1);
257 element->setNodeIndex(1, element->getNodeIndex(1) - 1);
258 element->setNodeIndex(2, element->getNodeIndex(2) - 1);
259 element->setNodeIndex(3, element->getNodeIndex(3) - 1);
260
261 element->setIndex(1); // default (no multibody meshing atm)
262 }
263
264 // Find parents for surface elements:
265 //------------------------------------
266 meshutils.findSurfaceElementParents(mesh);
267
268 // Find edges for surface elements:
269 //----------------------------------
270 meshutils.findSurfaceElementEdges(mesh);
271
272 // Compute normals for boundary elements:
273 //---------------------------------------
274 meshutils.findSurfaceElementNormals(mesh);
275
276 mesh->setDim(ngDim);
277 mesh->setCdim(ngDim);
278 }
279