src/escriptreader/MeshWithElements.cpp


/*******************************************************
*
* Copyright (c) 2003-2008 by University of Queensland
* Earth Systems Science Computational Center (ESSCC)
* http://www.uq.edu.au/esscc
*
* Primary Business: Queensland, Australia
* Licensed under the Open Software License version 3.0
* http://www.opensource.org/licenses/osl-3.0.php
*
*******************************************************/

//
// MeshWithElements.cpp
//
#include <escriptreader/MeshWithElements.h>
#include <escriptreader/ElementData.h>
#include <netcdf.hh>
#if HAVE_SILO
#include <silo.h>
#endif

using namespace std;

namespace EscriptReader {

//
//
//
MeshWithElements::MeshWithElements() : Mesh()
{
    name = "Elements";
}

//
//
//
MeshWithElements::MeshWithElements(const MeshWithElements& m) :
    Mesh(m)
{
    nodeTag = m.nodeTag;
    nodeGDOF = m.nodeGDOF;
    nodeGNI = m.nodeGNI;
    nodeGRDFI = m.nodeGRDFI;
    nodeGRNI = m.nodeGRNI;
    cells = new ElementData(*m.cells);
    faces = new ElementData(*m.faces);
    contacts = new ElementData(*m.contacts);
    points = new ElementData(*m.points);
}

//
//
//
MeshWithElements::~MeshWithElements()
{
    delete cells;
    delete faces;
    delete contacts;
    delete points;
}

//
// Returns a vector of strings containing Silo mesh names that have been
// written
//
StringVec MeshWithElements::getMeshNames() const
{
    StringVec res;
    res.push_back(name);
    StringVec tmpVec;
    tmpVec = cells->getMeshNames();
    res.insert(res.end(), tmpVec.begin(), tmpVec.end());
    tmpVec = faces->getMeshNames();
    res.insert(res.end(), tmpVec.begin(), tmpVec.end());
    tmpVec = contacts->getMeshNames();
    res.insert(res.end(), tmpVec.begin(), tmpVec.end());
    tmpVec = points->getMeshNames();
    res.insert(res.end(), tmpVec.begin(), tmpVec.end());
    return res;
}

Mesh* MeshWithElements::getMeshByName(const string name) const
{
    Mesh* ret = NULL;
    if (name == "Elements")
        ret = cells->fullMesh;
    else if (name == "ReducedElements")
        ret = cells->reducedMesh;
    else if (name == "FaceElements")
        ret = faces->fullMesh;
    else if (name == "ReducedFaceElements")
        ret = faces->reducedMesh;
    else if (name == "ContactElements")
        ret = contacts->fullMesh;
    else if (name == "ReducedContactElements")
        ret = contacts->reducedMesh;
    else if (name == "Points")
        ret = points->fullMesh;

    return ret;
}

//
// Returns a vector of strings containing Silo variable names that have
// been written
//
StringVec MeshWithElements::getVarNames() const
{
    StringVec res;
    res.push_back("Nodes_Id");
    res.push_back("Nodes_Tag");
    res.push_back("Nodes_gDOF");
    res.push_back("Nodes_gNI");
    res.push_back("Nodes_grDfI");
    res.push_back("Nodes_grNI");

    StringVec tmpVec;
    tmpVec = cells->getVarNames();
    res.insert(res.end(), tmpVec.begin(), tmpVec.end());
    tmpVec = faces->getVarNames();
    res.insert(res.end(), tmpVec.begin(), tmpVec.end());
    tmpVec = contacts->getVarNames();
    res.insert(res.end(), tmpVec.begin(), tmpVec.end());
    tmpVec = points->getVarNames();
    res.insert(res.end(), tmpVec.begin(), tmpVec.end());

    return res;
}

const IntVec& MeshWithElements::getVarDataByName(const string name) const
{
    if (name == "Nodes_Id")
        return nodeID;
    else if (name == "Nodes_Tag")
        return nodeTag;
    else if (name == "Nodes_gDOF")
        return nodeGDOF;
    else if (name == "Nodes_gNI")
        return nodeGNI;
    else if (name == "Nodes_grDfI")
        return nodeGRDFI;
    else if (name == "Nodes_grNI")
        return nodeGRNI;
    else if (name.compare(0, 9, "Elements_") == 0)
        return cells->getVarDataByName(name);
    else if (name.compare(0, 13, "FaceElements_") == 0)
        return faces->getVarDataByName(name);
    else if (name.compare(0, 16, "ContactElements_") == 0)
        return contacts->getVarDataByName(name);
    else if (name.compare(0, 7, "Points_") == 0)
        return points->getVarDataByName(name);
    else
        return *(IntVec*)(NULL);
}

ElementData* MeshWithElements::getElementsByName(const string name) const
{
    ElementData* ret = NULL;
    if (name == "Elements" || name == "ReducedElements")
        ret = cells;
    else if (name == "FaceElements" || name == "ReducedFaceElements")
        ret = faces;
    else if (name == "ContactElements" || name == "ReducedContactElements")
        ret = contacts;
    else if (name == "Points")
        ret = points;

    return ret;
}

//
// Reads mesh and element data from NetCDF file with given name
//
bool MeshWithElements::readFromNc(const string& filename)
{
    if (!Mesh::readFromNc(filename))
        return false;

    NcError ncerr(NcError::silent_nonfatal);
    NcFile* input;
    NcVar* var;
 
    input = new NcFile(filename.c_str());
    if (!input->is_valid()) {
        cerr << "Could not open input file " << filename << "." << endl;
        delete input;
        return false;
    }

    nodeTag.clear();
    nodeTag.insert(nodeTag.end(), numNodes, 0);
    var = input->get_var("Nodes_Tag");
    var->get(&nodeTag[0], numNodes);

    nodeGDOF.clear();
    nodeGDOF.insert(nodeGDOF.end(), numNodes, 0);
    var = input->get_var("Nodes_gDOF");
    var->get(&nodeGDOF[0], numNodes);

    nodeGNI.clear();
    nodeGNI.insert(nodeGNI.end(), numNodes, 0);
    var = input->get_var("Nodes_gNI");
    var->get(&nodeGNI[0], numNodes);

    nodeGRDFI.clear();
    nodeGRDFI.insert(nodeGRDFI.end(), numNodes, 0);
    var = input->get_var("Nodes_grDfI");
    var->get(&nodeGRDFI[0], numNodes);

    nodeGRNI.clear();
    nodeGRNI.insert(nodeGRNI.end(), numNodes, 0);
    var = input->get_var("Nodes_grNI");
    var->get(&nodeGRNI[0], numNodes);

    // Read all element types
    cells = new ElementData("Elements", this);
    cells->readFromNc(input);
    faces = new ElementData("FaceElements", this);
    faces->readFromNc(input);
    contacts = new ElementData("ContactElements", this);
    contacts->readFromNc(input);
    points = new ElementData("Points", this);
    points->readFromNc(input);

    delete input;
    return true;
}

//
//
//
void MeshWithElements::handleGhostZones(int ownIndex)
{
    cells->handleGhostZones(ownIndex);
    faces->handleGhostZones(ownIndex);
    contacts->handleGhostZones(ownIndex);
    points->handleGhostZones(ownIndex);
#ifdef _DEBUG
    cout << "block " << ownIndex << " has " << cells->getGhostCount()
         << " (" << cells->getReducedGhostCount() << ") ghost zones," << endl;
    cout << faces->getGhostCount() << " (" << faces->getReducedGhostCount()
         << ") ghost faces," << endl;
    cout << contacts->getGhostCount() << " ("
         << contacts->getReducedGhostCount() << ") ghost contacts," << endl;
    cout << points->getGhostCount() << " (" << points->getReducedGhostCount()
         << ") ghost points." << endl;
#endif
}

//
//
//
void MeshWithElements::removeGhostZones()
{
    cells->removeGhostZones();
    faces->removeGhostZones();
    contacts->removeGhostZones();
    points->removeGhostZones();
#ifdef _DEBUG
    cout << "After removing ghost zones there are" << endl;
    cout << "    " << numNodes << " Nodes, ";
    cout << cells->getCount() << " Elements, ";
    cout << faces->getCount() << " Face elements, ";
    cout << contacts->getCount() << " Contact elements, ";
    cout << points->getCount() << " Points left." << endl;
#endif
}

//
//
//
bool MeshWithElements::writeToSilo(DBfile* dbfile, const string& pathInSilo)
{
#if HAVE_SILO
    int ret;

    // Write meshes and zone-centered variables
    if (!cells->writeToSilo(dbfile, pathInSilo) ||
            !faces->writeToSilo(dbfile, pathInSilo) ||
            !contacts->writeToSilo(dbfile, pathInSilo) ||
            !points->writeToSilo(dbfile, pathInSilo))
        return false;

    if (!Mesh::writeToSilo(dbfile, pathInSilo))
        return false;
    
    if (pathInSilo != "") {
        ret = DBSetDir(dbfile, pathInSilo.c_str());
        if (ret != 0)
            return false;
    }

    string siloMeshName = getFullSiloName();

    // Write node-centered variables
    ret = DBPutUcdvar1(dbfile, "Nodes_Tag", siloMeshName.c_str(),
                (float*)&nodeTag[0], numNodes, NULL, 0, DB_INT,
                DB_NODECENT, NULL);
    if (ret == 0)
        ret = DBPutUcdvar1(dbfile, "Nodes_gDOF", siloMeshName.c_str(),
                (float*)&nodeGDOF[0], numNodes, NULL, 0, DB_INT,
                DB_NODECENT, NULL);
    if (ret == 0)
        ret = DBPutUcdvar1(dbfile, "Nodes_gNI", siloMeshName.c_str(),
                (float*)&nodeGNI[0], numNodes, NULL, 0, DB_INT,
                DB_NODECENT, NULL);
    if (ret == 0)
        ret = DBPutUcdvar1(dbfile, "Nodes_grDfI", siloMeshName.c_str(),
                (float*)&nodeGRDFI[0], numNodes, NULL, 0, DB_INT,
                DB_NODECENT, NULL);
    if (ret == 0)
        ret = DBPutUcdvar1(dbfile, "Nodes_grNI", siloMeshName.c_str(),
                (float*)&nodeGRNI[0], numNodes, NULL, 0, DB_INT,
                DB_NODECENT, NULL);

    DBSetDir(dbfile, "/");
    return (ret == 0);

#else // !HAVE_SILO
    return false;
#endif
}

} // namespace EscriptReader

1	caltinay	2183
2			/*******************************************************
3			*
4			* Copyright (c) 2003-2008 by University of Queensland
5			* Earth Systems Science Computational Center (ESSCC)
6			* http://www.uq.edu.au/esscc
7			*
8			* Primary Business: Queensland, Australia
9			* Licensed under the Open Software License version 3.0
10			* http://www.opensource.org/licenses/osl-3.0.php
11			*
12			*******************************************************/
13
14			//
15			// MeshWithElements.cpp
16			//
17			#include <escriptreader/MeshWithElements.h>
18			#include <escriptreader/ElementData.h>
19			#include <netcdf.hh>
20			#if HAVE_SILO
21			#include <silo.h>
22			#endif
23
24			using namespace std;
25
26	caltinay	2187	namespace EscriptReader {
27
28	caltinay	2183	//
29			//
30			//
31			MeshWithElements::MeshWithElements() : Mesh()
32			{
33			name = "Elements";
34			}
35
36			//
37			//
38			//
39			MeshWithElements::MeshWithElements(const MeshWithElements& m) :
40			Mesh(m)
41			{
42			nodeTag = m.nodeTag;
43			nodeGDOF = m.nodeGDOF;
44			nodeGNI = m.nodeGNI;
45			nodeGRDFI = m.nodeGRDFI;
46			nodeGRNI = m.nodeGRNI;
47			cells = new ElementData(*m.cells);
48			faces = new ElementData(*m.faces);
49			contacts = new ElementData(*m.contacts);
50			points = new ElementData(*m.points);
51			}
52
53			//
54			//
55			//
56			MeshWithElements::~MeshWithElements()
57			{
58			delete cells;
59			delete faces;
60			delete contacts;
61			delete points;
62			}
63
64			//
65			// Returns a vector of strings containing Silo mesh names that have been
66			// written
67			//
68			StringVec MeshWithElements::getMeshNames() const
69			{
70			StringVec res;
71			res.push_back(name);
72			StringVec tmpVec;
73			tmpVec = cells->getMeshNames();
74			res.insert(res.end(), tmpVec.begin(), tmpVec.end());
75			tmpVec = faces->getMeshNames();
76			res.insert(res.end(), tmpVec.begin(), tmpVec.end());
77			tmpVec = contacts->getMeshNames();
78			res.insert(res.end(), tmpVec.begin(), tmpVec.end());
79			tmpVec = points->getMeshNames();
80			res.insert(res.end(), tmpVec.begin(), tmpVec.end());
81			return res;
82			}
83
84			Mesh* MeshWithElements::getMeshByName(const string name) const
85			{
86			Mesh* ret = NULL;
87			if (name == "Elements")
88			ret = cells->fullMesh;
89			else if (name == "ReducedElements")
90			ret = cells->reducedMesh;
91			else if (name == "FaceElements")
92			ret = faces->fullMesh;
93			else if (name == "ReducedFaceElements")
94			ret = faces->reducedMesh;
95			else if (name == "ContactElements")
96			ret = contacts->fullMesh;
97			else if (name == "ReducedContactElements")
98			ret = contacts->reducedMesh;
99			else if (name == "Points")
100			ret = points->fullMesh;
101
102			return ret;
103			}
104
105			//
106			// Returns a vector of strings containing Silo variable names that have
107			// been written
108			//
109			StringVec MeshWithElements::getVarNames() const
110			{
111			StringVec res;
112			res.push_back("Nodes_Id");
113			res.push_back("Nodes_Tag");
114			res.push_back("Nodes_gDOF");
115			res.push_back("Nodes_gNI");
116			res.push_back("Nodes_grDfI");
117			res.push_back("Nodes_grNI");
118
119			StringVec tmpVec;
120			tmpVec = cells->getVarNames();
121			res.insert(res.end(), tmpVec.begin(), tmpVec.end());
122			tmpVec = faces->getVarNames();
123			res.insert(res.end(), tmpVec.begin(), tmpVec.end());
124			tmpVec = contacts->getVarNames();
125			res.insert(res.end(), tmpVec.begin(), tmpVec.end());
126			tmpVec = points->getVarNames();
127			res.insert(res.end(), tmpVec.begin(), tmpVec.end());
128
129			return res;
130			}
131
132			const IntVec& MeshWithElements::getVarDataByName(const string name) const
133			{
134			if (name == "Nodes_Id")
135			return nodeID;
136			else if (name == "Nodes_Tag")
137			return nodeTag;
138			else if (name == "Nodes_gDOF")
139			return nodeGDOF;
140			else if (name == "Nodes_gNI")
141			return nodeGNI;
142			else if (name == "Nodes_grDfI")
143			return nodeGRDFI;
144			else if (name == "Nodes_grNI")
145			return nodeGRNI;
146			else if (name.compare(0, 9, "Elements_") == 0)
147			return cells->getVarDataByName(name);
148			else if (name.compare(0, 13, "FaceElements_") == 0)
149			return faces->getVarDataByName(name);
150			else if (name.compare(0, 16, "ContactElements_") == 0)
151			return contacts->getVarDataByName(name);
152			else if (name.compare(0, 7, "Points_") == 0)
153			return points->getVarDataByName(name);
154			else
155			return (IntVec)(NULL);
156			}
157
158			ElementData* MeshWithElements::getElementsByName(const string name) const
159			{
160			ElementData* ret = NULL;
161			if (name == "Elements" \|\| name == "ReducedElements")
162			ret = cells;
163			else if (name == "FaceElements" \|\| name == "ReducedFaceElements")
164			ret = faces;
165			else if (name == "ContactElements" \|\| name == "ReducedContactElements")
166			ret = contacts;
167			else if (name == "Points")
168			ret = points;
169
170			return ret;
171			}
172
173			//
174			// Reads mesh and element data from NetCDF file with given name
175			//
176			bool MeshWithElements::readFromNc(const string& filename)
177			{
178			if (!Mesh::readFromNc(filename))
179			return false;
180
181			NcError ncerr(NcError::silent_nonfatal);
182			NcFile* input;
183			NcVar* var;
184
185			input = new NcFile(filename.c_str());
186			if (!input->is_valid()) {
187			cerr << "Could not open input file " << filename << "." << endl;
188			delete input;
189			return false;
190			}
191
192			nodeTag.clear();
193			nodeTag.insert(nodeTag.end(), numNodes, 0);
194			var = input->get_var("Nodes_Tag");
195			var->get(&nodeTag[0], numNodes);
196
197			nodeGDOF.clear();
198			nodeGDOF.insert(nodeGDOF.end(), numNodes, 0);
199			var = input->get_var("Nodes_gDOF");
200			var->get(&nodeGDOF[0], numNodes);
201
202			nodeGNI.clear();
203			nodeGNI.insert(nodeGNI.end(), numNodes, 0);
204			var = input->get_var("Nodes_gNI");
205			var->get(&nodeGNI[0], numNodes);
206
207			nodeGRDFI.clear();
208			nodeGRDFI.insert(nodeGRDFI.end(), numNodes, 0);
209			var = input->get_var("Nodes_grDfI");
210			var->get(&nodeGRDFI[0], numNodes);
211
212			nodeGRNI.clear();
213			nodeGRNI.insert(nodeGRNI.end(), numNodes, 0);
214			var = input->get_var("Nodes_grNI");
215			var->get(&nodeGRNI[0], numNodes);
216
217			// Read all element types
218			cells = new ElementData("Elements", this);
219			cells->readFromNc(input);
220			faces = new ElementData("FaceElements", this);
221			faces->readFromNc(input);
222			contacts = new ElementData("ContactElements", this);
223			contacts->readFromNc(input);
224			points = new ElementData("Points", this);
225			points->readFromNc(input);
226
227			delete input;
228			return true;
229			}
230
231			//
232			//
233			//
234			void MeshWithElements::handleGhostZones(int ownIndex)
235			{
236			cells->handleGhostZones(ownIndex);
237			faces->handleGhostZones(ownIndex);
238			contacts->handleGhostZones(ownIndex);
239			points->handleGhostZones(ownIndex);
240			#ifdef _DEBUG
241			cout << "block " << ownIndex << " has " << cells->getGhostCount()
242			<< " (" << cells->getReducedGhostCount() << ") ghost zones," << endl;
243			cout << faces->getGhostCount() << " (" << faces->getReducedGhostCount()
244			<< ") ghost faces," << endl;
245			cout << contacts->getGhostCount() << " ("
246			<< contacts->getReducedGhostCount() << ") ghost contacts," << endl;
247			cout << points->getGhostCount() << " (" << points->getReducedGhostCount()
248			<< ") ghost points." << endl;
249			#endif
250			}
251
252			//
253			//
254			//
255			void MeshWithElements::removeGhostZones()
256			{
257			cells->removeGhostZones();
258			faces->removeGhostZones();
259			contacts->removeGhostZones();
260			points->removeGhostZones();
261			#ifdef _DEBUG
262			cout << "After removing ghost zones there are" << endl;
263			cout << " " << numNodes << " Nodes, ";
264			cout << cells->getCount() << " Elements, ";
265			cout << faces->getCount() << " Face elements, ";
266			cout << contacts->getCount() << " Contact elements, ";
267			cout << points->getCount() << " Points left." << endl;
268			#endif
269			}
270
271			//
272			//
273			//
274			bool MeshWithElements::writeToSilo(DBfile* dbfile, const string& pathInSilo)
275			{
276			#if HAVE_SILO
277			int ret;
278
279			// Write meshes and zone-centered variables
280			if (!cells->writeToSilo(dbfile, pathInSilo) \|\|
281			!faces->writeToSilo(dbfile, pathInSilo) \|\|
282			!contacts->writeToSilo(dbfile, pathInSilo) \|\|
283			!points->writeToSilo(dbfile, pathInSilo))
284			return false;
285
286			if (!Mesh::writeToSilo(dbfile, pathInSilo))
287			return false;
288
289			if (pathInSilo != "") {
290			ret = DBSetDir(dbfile, pathInSilo.c_str());
291			if (ret != 0)
292			return false;
293			}
294
295			string siloMeshName = getFullSiloName();
296
297			// Write node-centered variables
298			ret = DBPutUcdvar1(dbfile, "Nodes_Tag", siloMeshName.c_str(),
299			(float*)&nodeTag[0], numNodes, NULL, 0, DB_INT,
300			DB_NODECENT, NULL);
301			if (ret == 0)
302			ret = DBPutUcdvar1(dbfile, "Nodes_gDOF", siloMeshName.c_str(),
303			(float*)&nodeGDOF[0], numNodes, NULL, 0, DB_INT,
304			DB_NODECENT, NULL);
305			if (ret == 0)
306			ret = DBPutUcdvar1(dbfile, "Nodes_gNI", siloMeshName.c_str(),
307			(float*)&nodeGNI[0], numNodes, NULL, 0, DB_INT,
308			DB_NODECENT, NULL);
309			if (ret == 0)
310			ret = DBPutUcdvar1(dbfile, "Nodes_grDfI", siloMeshName.c_str(),
311			(float*)&nodeGRDFI[0], numNodes, NULL, 0, DB_INT,
312			DB_NODECENT, NULL);
313			if (ret == 0)
314			ret = DBPutUcdvar1(dbfile, "Nodes_grNI", siloMeshName.c_str(),
315			(float*)&nodeGRNI[0], numNodes, NULL, 0, DB_INT,
316			DB_NODECENT, NULL);
317
318			DBSetDir(dbfile, "/");
319			return (ret == 0);
320
321			#else // !HAVE_SILO
322			return false;
323			#endif
324			}
325
326	caltinay	2187	} // namespace EscriptReader
327