00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #include <ScalarFunctionReaderVTK.hpp>
00021 #include <ScalarFunctionBase.hpp>
00022
00023 #include <ErrorHandler.hpp>
00024 #include <Mesh.hpp>
00025
00026 #include <Structured3DMesh.hpp>
00027
00028 #include <ScalarDiscretizationTypeBase.hpp>
00029 #include <Vector.hpp>
00030 #include <FEMFunctionBase.hpp>
00031
00032 #include <FEMFunctionBuilder.hpp>
00033
00034 #include <XMLFileReader.hpp>
00035 #include <XMLTree.hpp>
00036
00037 #include <XMLContentPosition.hpp>
00038
00039 #include <EndianConverter.hpp>
00040
00041 #include <fstream>
00042 #include <sstream>
00043
00044 ConstReferenceCounting<ScalarFunctionBase>
00045 ScalarFunctionReaderVTK::
00046 getFunction() const
00047 {
00048 ffout(3) << "Reading file '" << __filename << "'\n";
00049 std::ifstream fin(__filename.c_str());
00050
00051 if (not(fin)) {
00052 throw ErrorHandler(__FILE__,__LINE__,
00053 "error: cannot open file '"+__filename+"'",
00054 ErrorHandler::normal);
00055 }
00056
00057 XMLFileReader::instance().read(__filename);
00058
00059 const XMLTree& xmlTree = *XMLFileReader::instance().xmlTree();
00060 xmlTree.check();
00061
00062
00063 ReferenceCounting<XMLTag> vtkfile = xmlTree.findTag(">VTKFile");
00064 ReferenceCounting<XMLAttribute> filetype = vtkfile->findAttribute("type");
00065 ffout(3) << " - file type is '" << filetype->value() << "'\n";
00066
00067 Vector<real_t> v;
00068 size_t numberOfCells = 0;
00069 size_t numberOfPoints = 0;
00070 size_t numberOfComponents = 0;
00071
00072 if (filetype->value() == "UnstructuredGrid") {
00073 ReferenceCounting<XMLTag> piece = xmlTree.findTag(">VTKFile>UnstructuredGrid>Piece");
00074
00075
00076 {
00077 ReferenceCounting<XMLAttribute> attribute = piece->findAttribute("NumberOfCells");
00078 std::istringstream is(attribute->value());
00079 is >> numberOfCells;
00080 }
00081
00082 {
00083 ReferenceCounting<XMLAttribute> attribute = piece->findAttribute("NumberOfPoints");
00084 std::istringstream is(attribute->value());
00085 is >> numberOfPoints;
00086 }
00087 } else if (filetype->value() == "ImageData") {
00088 ReferenceCounting<XMLTag> piece = xmlTree.findTag(">VTKFile>ImageData>Piece");
00089
00090
00091 {
00092 ReferenceCounting<XMLAttribute> attribute = piece->findAttribute("Extent");
00093 std::istringstream is(attribute->value());
00094 TinyVector<3,size_t> extent0;
00095 TinyVector<3,size_t> extent1;
00096 is >> extent0[0] >> extent1[0] >> extent0[1] >> extent1[1] >> extent0[2] >> extent1[2];
00097 numberOfCells = (extent1[0]-extent0[0])*(extent1[1]-extent0[1])*(extent1[2]-extent0[2]);
00098 numberOfPoints = (1+extent1[0]-extent0[0])*(1+extent1[1]-extent0[1])*(1+extent1[2]-extent0[2]);
00099 }
00100 } else {
00101 throw ErrorHandler(__FILE__,__LINE__,
00102 "cannot read filetype '"+filetype->value()+"'",
00103 ErrorHandler::unexpected);
00104 }
00105
00106
00107 if ((numberOfPoints != __mesh->numberOfVertices()) or
00108 (numberOfCells != __mesh->numberOfCells())) {
00109 throw ErrorHandler(__FILE__,__LINE__,
00110 "mesh in file '"+__filename+"' is different from 'read' argument",
00111 ErrorHandler::normal);
00112 }
00113
00114 ReferenceCounting<XMLTag> dataArray = 0;
00115
00116
00117 if (xmlTree.hasTag(">VTKFile>"+filetype->value()+">Piece>PointData>DataArray")) {
00118 XMLTree::const_range pointDataRange = xmlTree.findTagRange(">VTKFile>"+filetype->value()+">Piece>PointData>DataArray");
00119 for (XMLTree::TagList::const_iterator iDataArray = pointDataRange.first;
00120 iDataArray != pointDataRange.second; ++iDataArray) {
00121 ReferenceCounting<XMLAttribute> name = iDataArray->second->findAttribute("Name");
00122 if (name->value() == __functionName) {
00123 dataArray = iDataArray->second;
00124 ReferenceCounting<XMLAttribute> components = dataArray->findAttribute("NumberOfComponents","1");
00125 std::istringstream is(components->value());
00126 is >> numberOfComponents;
00127 break;
00128 }
00129 }
00130 }
00131
00132 if (dataArray == 0) {
00133 if (xmlTree.hasTag(">VTKFile>"+filetype->value()+">Piece>CellData>DataArray")) {
00134 XMLTree::const_range cellDataRange = xmlTree.findTagRange(">VTKFile>"+filetype->value()+">Piece>CellData>DataArray");
00135 for (XMLTree::TagList::const_iterator iDataArray = cellDataRange.first;
00136 iDataArray != cellDataRange.second; ++iDataArray) {
00137 ReferenceCounting<XMLAttribute> name = iDataArray->second->findAttribute("Name");
00138 if (name->value() == __functionName) {
00139 dataArray = iDataArray->second;
00140 ReferenceCounting<XMLAttribute> components = dataArray->findAttribute("NumberOfComponents","1");
00141 std::istringstream is(components->value());
00142 is >> numberOfComponents;
00143 break;
00144 }
00145 }
00146 }
00147 }
00148
00149 if (dataArray == 0) {
00150 throw ErrorHandler(__FILE__,__LINE__,
00151 "cannot find field '"+__functionName+"' in '"+__filename+"'",
00152 ErrorHandler::normal);
00153 }
00154
00155 ConstReferenceCounting<XMLAttribute> format = dataArray->findAttribute("format");
00156 ffout(3) << " - format for \'" << __functionName << "\' is: " << format->value() << '\n';
00157
00158 if (__componentNumber>=numberOfComponents) {
00159 throw ErrorHandler(__FILE__,__LINE__,
00160 "The function '"+__functionName+"' of the file '"+__filename
00161 +"' does only contain "+stringify(__componentNumber)+" components.\n"
00162 +"Cannot read component number "+stringify(__componentNumber),
00163 ErrorHandler::normal);
00164 }
00165
00166 v.resize(numberOfPoints);
00167
00168 if (format->value() == "appended") {
00169 ConstReferenceCounting<XMLAttribute> offset = dataArray->findAttribute("offset");
00170 if (offset == 0) {
00171 throw ErrorHandler(__FILE__,__LINE__,
00172 "cannot find offset for '"+__functionName+"' in '"+__filename+"'",
00173 ErrorHandler::normal);
00174 }
00175
00176 ReferenceCounting<XMLTag> appendedData = xmlTree.findTag(">VTKFile>AppendedData");
00177
00178 if (appendedData == 0) {
00179 throw ErrorHandler(__FILE__,__LINE__,
00180 "cannot find AppendedData to read '"+__functionName+"' in '"+__filename+"'",
00181 ErrorHandler::normal);
00182 }
00183 ConstReferenceCounting<XMLContentBase> content = appendedData->content();
00184 if (content == 0) {
00185 throw ErrorHandler(__FILE__,__LINE__,
00186 "cannot find AppendedData content for '"+__functionName+"' in '"+__filename+"'",
00187 ErrorHandler::normal);
00188 }
00189
00190 ReferenceCounting<XMLAttribute> endianness = vtkfile->findAttribute("byte_order");
00191 ffout(3) << " - endianness is '" << endianness->value() << "'\n";
00192 bool littleEndian = true;
00193 if (endianness->value() == "LittleEndian") {
00194 littleEndian = true;
00195 } else if (endianness->value() == "BigEndian") {
00196 littleEndian = false;
00197 } else {
00198 throw ErrorHandler(__FILE__,__LINE__,
00199 "unknown endianness type",
00200 ErrorHandler::normal);
00201 }
00202
00203 switch(content->type()) {
00204 case XMLContentBase::filePosition: {
00205 const XMLContentPosition& position = dynamic_cast<const XMLContentPosition&>(*content);
00206
00207 fin.seekg(position.position(),std::ios_base::beg);
00208 fin.unget();
00209
00210 char c =0;
00211 fin.get(c);
00212 if(c != '_') {
00213 throw ErrorHandler(__FILE__,__LINE__,
00214 "'_' character disappeared since parser was called!",
00215 ErrorHandler::unexpected);
00216 }
00217
00218 std::stringstream s;
00219 s << offset->value() << std::ends;
00220 int offsetValue = 0;
00221 s >> offsetValue;
00222
00223 fin.seekg(offsetValue,std::ios_base::cur);
00224
00225 int dataSize = 0;
00226 fin.read(reinterpret_cast<char*>(&dataSize),sizeof(int));
00227 if (littleEndian) {
00228 fromLittleEndian(dataSize);
00229 } else {
00230 fromBigEndian(dataSize);
00231 }
00232
00233 if (static_cast<size_t>(dataSize) != sizeof(real_t)*numberOfPoints*numberOfComponents) {
00234 throw ErrorHandler(__FILE__,__LINE__,
00235 "data storage inconsistency while reading function '"+__functionName+"' of file '"+__filename+"'",
00236 ErrorHandler::normal);
00237 }
00238
00239 Vector<real_t> allValues(numberOfComponents*numberOfPoints);
00240 fin.read(reinterpret_cast<char*>(&(allValues[0])),sizeof(real_t)*allValues.size());
00241 if (littleEndian) {
00242 fromLittleEndian(allValues);
00243 } else {
00244 fromBigEndian(allValues);
00245 }
00246
00247 switch (__mesh->type()) {
00248 case Mesh::cartesianHexahedraMesh: {
00249 const Structured3DMesh& mesh = dynamic_cast<const Structured3DMesh&>(*__mesh);
00250
00251 size_t number=0;
00252 for (size_t k=0; k<mesh.shape().nz(); ++k) {
00253 for (size_t j=0; j<mesh.shape().ny(); ++j) {
00254 for (size_t i=0; i<mesh.shape().nx(); ++i) {
00255 v[mesh.shape()(i,j,k)] = allValues[number*numberOfComponents+__componentNumber];
00256 number++;
00257 }
00258 }
00259 }
00260 break;
00261 }
00262 default: {
00263 for (size_t i=0; i<numberOfPoints; ++i) {
00264 v[i] = allValues[numberOfComponents*i+__componentNumber];
00265 }
00266 }
00267 }
00268 break;
00269 }
00270 default: {
00271 throw ErrorHandler(__FILE__,__LINE__,
00272 "Bad content for '"+__functionName+"' in '"+__filename+"'",
00273 ErrorHandler::normal);
00274
00275 }
00276 }
00277
00278 } else if (format->value() == "ascii") {
00279 ConstReferenceCounting<XMLContentBase> content = dataArray->content();
00280
00281 if (content == 0) {
00282 throw ErrorHandler(__FILE__,__LINE__,
00283 "cannot find content for '"+__functionName+"' in '"+__filename+"'",
00284 ErrorHandler::normal);
00285 }
00286
00287 switch(content->type()) {
00288 case XMLContentBase::filePosition: {
00289 const XMLContentPosition& position = dynamic_cast<const XMLContentPosition&>(*content);
00290
00291 fin.seekg(position.position(),std::ios_base::beg);
00292 fin.unget();
00293 switch (__mesh->type()) {
00294 case Mesh::cartesianHexahedraMesh: {
00295 const Structured3DMesh& mesh = dynamic_cast<const Structured3DMesh&>(*__mesh);
00296 Vector<real_t> componentsValue(numberOfComponents);
00297
00298 for (size_t k=0; k<mesh.shape().nz(); ++k) {
00299 for (size_t j=0; j<mesh.shape().ny(); ++j) {
00300 for (size_t i=0; i<mesh.shape().nx(); ++i) {
00301 for (size_t l=0; l<componentsValue.size(); ++l) {
00302 fin >> componentsValue[l];
00303 }
00304 v[mesh.shape()(i,j,k)] = componentsValue[__componentNumber];
00305 }
00306 }
00307 }
00308 break;
00309 }
00310 default: {
00311 Vector<real_t> componentsValue(numberOfComponents);
00312 for (size_t i=0; i<numberOfPoints; ++i) {
00313 for (size_t j=0; j<componentsValue.size(); ++j) {
00314 fin >> componentsValue[j];
00315 }
00316 v[i] = componentsValue[__componentNumber];
00317 }
00318 }
00319 }
00320 break;
00321 }
00322 default: {
00323 throw ErrorHandler(__FILE__,__LINE__,
00324 "Bad content for '"+__functionName+"' in '"+__filename+"'",
00325 ErrorHandler::normal);
00326
00327 }
00328 }
00329 } else {
00330 throw ErrorHandler(__FILE__,__LINE__,
00331 "Format '"+format->value()+"' is not supported",
00332 ErrorHandler::normal);
00333 }
00334
00335 FEMFunctionBuilder builder;
00336 ScalarDiscretizationTypeFEM d(ScalarDiscretizationTypeBase::lagrangianFEM1);
00337 builder.build(d, __mesh, v);
00338 return builder.getBuiltScalarFunction();
00339 }
00340
00341 ScalarFunctionReaderVTK::
00342 ScalarFunctionReaderVTK(const std::string& filename,
00343 ConstReferenceCounting<Mesh> mesh,
00344 const std::string& functionName,
00345 const int& componentNumber)
00346 : ScalarFunctionReaderBase(filename, mesh, functionName, componentNumber)
00347 {
00348 ;
00349 }
00350
00351 ScalarFunctionReaderVTK::
00352 ScalarFunctionReaderVTK(const ScalarFunctionReaderVTK& f)
00353 : ScalarFunctionReaderBase(f)
00354 {
00355 ;
00356 }
00357
00358 ScalarFunctionReaderVTK::
00359 ~ScalarFunctionReaderVTK()
00360 {
00361 ;
00362 }
