00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #ifndef MECHSYS_FEM_OUTPUT_H
00023 #define MECHSYS_FEM_OUTPUT_H
00024
00025 #ifdef HAVE_CONFIG_H
00026 #include "config.h"
00027 #else
00028 #ifndef REAL
00029 #define REAL double
00030 #endif
00031 #endif
00032
00033 #define OUTPUT_MAX_NODES 20
00034 #define OUTPUT_MAX_ELEMS 10
00035
00036 #include <iostream>
00037 #include <sstream>
00038 #include <fstream>
00039 #include <string>
00040 #include <map>
00041
00042 #include "util/exception.h"
00043 #include "util/numstreams.h"
00044 #include "fem/inputdata.h"
00045 #include "fem/data.h"
00046 #include "fem/solver/intsolverdata.h"
00047 #include "linalg/matrix.h"
00048
00049 using Util::_4;
00050
00051 namespace FEM
00052 {
00053
00054 class Output
00055 {
00056 friend std::ostream & operator<< (std::ostream & os, FEM::Output const & OU);
00057 public:
00058 Output() : _n_out_nodes(0), _n_out_elems(0), _notify(NULL), _pt2object(NULL), _pt2wrapper(NULL) {}
00059 ~Output();
00060 void Initialize (FEM::InputData const * idata, FEM::Data const * data);
00061 void OutputIncrement (int iStage, int increment);
00062 void OutputIncrement (int iStage, int increment, IntSolverData const & ISD);
00063 void OutputStage (int iStage);
00064 void SetNotifyCallBack (void (*pFun)(int, int, IntSolverData const &)) { _notify=pFun; }
00065 void SetNotifyCallBack (void* pt2Object, void (*pt2Wrapper)(void* pt2Object, int, int) )
00066 {
00067 _pt2object = pt2Object;
00068 _pt2wrapper = pt2Wrapper;
00069 }
00070 void OpenFiles();
00071 void CloseFiles();
00072 private:
00073 FEM::InputData const * _idata;
00074 FEM::Data const * _data;
00075 int _n_out_nodes;
00076 int _n_out_elems;
00077 std::ofstream _a_node_files[OUTPUT_MAX_NODES];
00078 std::ofstream _a_elem_files[OUTPUT_MAX_ELEMS];
00079
00080 void VTKCreateFile(String const & Filename);
00081 void VTKAppendData(String const & Filename, int iStage, int iInc);
00082 void _write_stage_opendx(size_t iStage, Array<String> & Labels, LinAlg::Matrix<REAL> & DataValues);
00083 void _write_stage_vtk (size_t iStage, Array<String> & Labels, LinAlg::Matrix<REAL> & DataValues);
00084 private:
00085 void (*_notify)(int iStage, int increment, IntSolverData const & ISD);
00086 void * _pt2object;
00087 void (*_pt2wrapper)(void* pt2Object, int iStage, int increment);
00088
00089 };
00090
00091 std::ostream & operator<< (std::ostream & os, FEM::Output const & OU)
00092 {
00093 os << (*OU._data) << std::endl;
00094 return os;
00095 }
00096
00097
00099
00100
00101 inline void Output::Initialize(FEM::InputData const * idata, FEM::Data const * data)
00102 {
00103 _idata = idata;
00104 _data = data;
00105 _n_out_nodes = _idata->outNODES.size();
00106 _n_out_elems = _idata->outELEMS.size();
00107
00108
00109 if (_n_out_nodes>OUTPUT_MAX_NODES)
00110 throw new Fatal(_("Output::Initialize: Number of nodes for output too big <max=%d>"),OUTPUT_MAX_NODES);
00111
00112
00113 for (int i_node=0; i_node<_n_out_nodes; ++i_node)
00114 {
00115 int id=_idata->outNODES[i_node]-1;
00116 if (!(id>=0 && id<static_cast<int>(data->Nodes.size())))
00117 throw new Fatal(_("Output::Initialize: Number of node < %d > to output is out of range <max=%d>"), id+1, data->Nodes.size());
00118 }
00119
00120
00121 if (_n_out_elems>OUTPUT_MAX_ELEMS)
00122 throw new Fatal(_("Output::Initialize: Number of elements for output too big <max=%d>"),OUTPUT_MAX_ELEMS);
00123
00124
00125 for (int i_elem=0; i_elem<_n_out_elems; ++i_elem)
00126 {
00127 int id=_idata->outELEMS[i_elem]-1;
00128 if (!(id>=0 && id<static_cast<int>(data->Elements.size())))
00129 throw new Fatal(_("Output::Initialize: Number of element < %d > to output is out of range <max=%d>"), id+1, data->Elements.size());
00130 }
00131
00132
00133 }
00134
00135 Output::~Output()
00136 {
00137 CloseFiles();
00138 }
00139
00140 void Output::OutputIncrement(int iStage, int increment)
00141 {
00142 IntSolverData isd;
00143 OutputIncrement(iStage, increment, isd);
00144 }
00145
00146 void Output::OutputIncrement(int iStage, int increment, IntSolverData const & ISD)
00147 {
00148
00149 if (increment==-1) VTKCreateFile(_idata->fnOutVTK);
00150 else
00151 {
00152 if (_idata->vtkOutAllIncs) VTKAppendData(_idata->fnOutVTK, iStage, increment);
00153 else
00154 if (increment==_idata->nDIV[iStage]-1)
00155 VTKAppendData(_idata->fnOutVTK, iStage, increment);
00156 }
00157
00158
00159 for (int i_node=0; i_node<_n_out_nodes; ++i_node)
00160 {
00161 if (_data->Nodes[_idata->outNODES[i_node]-1].Refs()>0)
00162 {
00163 if (increment==-1)
00164 _a_node_files[i_node] << _6()<< "Inc" << _data->Nodes[ _idata->outNODES[i_node]-1 ].Out(true);
00165 _a_node_files[i_node] << _6()<< increment+1 << _data->Nodes[ _idata->outNODES[i_node]-1 ].Out();
00166 }
00167 }
00168
00169
00170 for (int i_elem=0; i_elem<_n_out_elems; ++i_elem)
00171 {
00172 if (_data->Elements[_idata->outELEMS[i_elem]-1]->IsActive())
00173 {
00174 if (increment==-1)
00175 _a_elem_files[i_elem] << _6()<< "Inc" << _data->Elements[ _idata->outELEMS[i_elem]-1 ]->OutCenter(true);
00176 _a_elem_files[i_elem] << _6()<< increment+1 << _data->Elements[ _idata->outELEMS[i_elem]-1 ]->OutCenter();
00177 }
00178 }
00179
00180 if (_notify !=NULL) (*_notify)(iStage, increment, ISD);
00181 if (_pt2object!=NULL) { if (_pt2wrapper!=NULL) (*_pt2wrapper)(_pt2object, iStage, increment); }
00182
00183 }
00184
00185 void Output::OutputStage(int iStage)
00186 {
00187 if (_idata->DoOutSTAGE[iStage]==0) return;
00188
00189 int n_nodes = _data->Nodes.size();
00190 Array<Element *> const & elements = _data->Elements;
00191 Array<Element *> active_elements;
00192 int n_elements = elements.size();
00193 int n_active_elems = 0;
00194
00195
00196
00197
00198 for (int i_elem=0; i_elem<n_elements; ++i_elem)
00199 {
00200 String ele_name = elements[i_elem]->Name();
00201 if (elements[i_elem]->IsActive() == false) continue;
00202 if (ele_name.substr(0,3) == "Lin") continue;
00203 active_elements.push_back(elements[i_elem]);
00204 n_active_elems++;
00205 }
00206
00207 int const MAX_DATA_COMP = 20;
00208 LinAlg::Matrix<REAL > values(n_nodes, MAX_DATA_COMP);
00209 LinAlg::Matrix<size_t> refs (n_nodes, MAX_DATA_COMP);
00210 values.SetValues(0.0);
00211 values.SetValues(0 );
00212
00213 std::map<String, int> index_map;
00214 Array<String> all_labels;
00215 int total_comps = 0;
00216
00217
00218 for (int i_elem=0; i_elem<n_active_elems; ++i_elem)
00219 {
00220 LinAlg::Matrix<REAL> elem_values;
00221 Array<String> elem_labels;
00222 active_elements[i_elem]->OutNodes(elem_values, elem_labels);
00223 int n_labels = elem_labels.size();
00224 Array<Node*> const & connects = active_elements[i_elem]->Connects();
00225 int n_elem_nodes = active_elements[i_elem]->nNodes();
00226 for (int j_label = 0; j_label < n_labels; ++j_label)
00227 {
00228 String & current_label = elem_labels[j_label];
00229 if (index_map.find(elem_labels[j_label]) == index_map.end())
00230 {
00231 total_comps++;
00232 all_labels.push_back(elem_labels[j_label]);
00233 index_map[current_label] = total_comps-1;
00234 }
00235
00236 int index = index_map[current_label];
00237
00238 for (int j_node=0; j_node<n_elem_nodes; ++j_node)
00239 {
00240 int node_number = connects[j_node]->Number();
00241 values(node_number,index) += elem_values(j_node, j_label);
00242 refs (node_number,index) ++;
00243 }
00244 }
00245 }
00246 assert(total_comps<MAX_DATA_COMP);
00247
00248
00249 for (int i=0; i<n_nodes; i++)
00250 for (int j=0; j<total_comps; j++)
00251 {
00252 if (refs(i,j) != 0)
00253 values(i,j) /= refs(i,j);
00254 else
00255 values(i,j) = 0.0;
00256 }
00257
00258 _write_stage_vtk (iStage, all_labels, values);
00259 _write_stage_opendx(iStage, all_labels, values);
00260
00261 }
00262
00263 inline void Output::_write_stage_opendx(size_t iStage, Array<String> & Labels, LinAlg::Matrix<REAL> & DataValues)
00264 {
00265
00266 std::ofstream ofile;
00267 std::ostringstream oss;
00268 oss << _idata->fnOutSTAGES[iStage].GetSTL().c_str() << ".dx";
00269 ofile.open(oss.str().c_str(), std::ios::out);
00270
00271
00272 Array<Node> const & N = _data->Nodes;
00273 Array<Element*> const & All_E = _data->Elements;
00274 Array<Element*> E;
00275
00276
00277 size_t n_nodes = N.size();
00278 size_t n_all_elems = All_E.size();
00279 size_t n_elems = 0;
00280 size_t n_comps = Labels.size();
00281
00282
00283
00284
00285 for (size_t i=0; i<n_all_elems; ++i)
00286 {
00287 String ele_name =All_E[i]->Name();
00288 if (All_E[i]->IsActive() == false) continue;
00289 if (ele_name.substr(0,3) == "Lin") continue;
00290 E.push_back(All_E[i]);
00291 n_elems++;
00292 }
00293
00294
00295
00296
00297
00298 ofile << "object \"nodes\" class array type float rank 1 shape 3 items " << n_nodes << " data follows" << std::endl;
00299 for (size_t i=0; i<n_nodes; ++i) { ofile << _8s()<<N[i].X() << _8s() << N[i].Y() << _8s() << N[i].Z() << std::endl; }
00300 ofile << std::endl;
00301 ofile << "object \"elements\" class array type int rank 1 shape ";
00302 if (E[0]->nNodes() != 20) ofile << E[0]->nNodes();
00303 else ofile << 8;
00304 ofile << " items " << n_elems << " data follows" << std::endl;
00305 for (size_t i=0; i<n_elems; i++) {
00306 size_t n_nodes = E[i]->nNodes();
00307 Array<Node*> const & connects = \
00308 E[i]->Connects();
00309 if (n_nodes==8) {
00310 ofile << _6() << connects[0]->Number();
00311 ofile << _6() << connects[1]->Number();
00312 ofile << _6() << connects[3]->Number();
00313 ofile << _6() << connects[2]->Number();
00314 ofile << _6() << connects[4]->Number();
00315 ofile << _6() << connects[5]->Number();
00316 ofile << _6() << connects[7]->Number();
00317 ofile << _6() << connects[6]->Number(); }
00318 else if (n_nodes==20) {
00319 ofile << _6() << connects[0] ->Number();
00320 ofile << _6() << connects[2] ->Number();
00321 ofile << _6() << connects[6] ->Number();
00322 ofile << _6() << connects[4] ->Number();
00323 ofile << _6() << connects[12]->Number();
00324 ofile << _6() << connects[14]->Number();
00325 ofile << _6() << connects[18]->Number();
00326 ofile << _6() << connects[16]->Number(); }
00327 else {
00328 for (int j=0; j<E[i]->nNodes(); ++j) { ofile << E[i]->Connects()[j]->Number() << " "; } }
00329 ofile << std::endl; }
00330 ofile << "attribute \"element type\" string \"cubes\"" << std::endl;
00331 ofile << "attribute \"ref\" string \"positions\"" << std::endl;
00332 ofile << std::endl;
00333 for (size_t i=0; i<n_comps; i++) { ofile << "object \"" << Labels[i] << "\" class array type float rank 0 items " << n_nodes << " data follows" << std::endl;
00334 for (size_t j=0; j<n_nodes; ++j) { ofile << _8s() << DataValues(j,i) << std::endl; } }
00335 ofile << std::endl;
00336 ofile << "object \"_field\" class field" << std::endl;
00337 ofile << "component \"positions\" value \"nodes\"" << std::endl;
00338 ofile << "component \"connections\" value \"elements\"" << std::endl;
00339 ofile << "component \"data\" value \"" << Labels[0] << "\"" << std::endl;
00340
00341 ofile.close();
00342
00343 }
00344
00345 inline void Output::_write_stage_vtk(size_t iStage, Array<String> & Labels, LinAlg::Matrix<REAL> & DataValues)
00346 {
00347
00348 std::ofstream ofile;
00349 std::ostringstream oss;
00350 oss << _idata->fnOutSTAGES[iStage].GetSTL().c_str() << ".vtk";
00351 ofile.open(oss.str().c_str(), std::ios::out);
00352
00353
00354 Array<Node> const & N = _data->Nodes;
00355 Array<Element*> const & E = _data->Elements;
00356
00357
00358 size_t n_nodes = N.size();
00359 size_t n_all_elems = E.size();
00360 size_t n_elems = 0;
00361 size_t n_comps = Labels.size();
00362
00363
00364 Array<Element * > active_E;
00365 for (size_t i=0; i <n_all_elems; ++i)
00366 {
00367 String ele_name = E[i]->Name();
00368 if (E[i]->IsActive() == false) continue;
00369 if (ele_name.substr(0,3) == "Lin") continue;
00370 active_E.push_back(E[i]);
00371 n_elems++;
00372 }
00373
00374
00375 int n_data = 0;
00376 for (size_t i=0; i<n_elems; ++i) n_data += 1 + active_E[i]->nNodes();
00377
00378
00379 const String DUX = "Dux";
00380 const String DUY = "Duy";
00381 const String DUZ = "Duz";
00382
00383
00384 ofile << "# vtk DataFile Version 3.0" << std::endl;
00385 ofile << "MechSys/FEM - " << _idata->projNAME << std::endl;
00386 ofile << "ASCII" << std::endl;
00387 ofile << "DATASET UNSTRUCTURED_GRID" << std::endl;
00388 ofile << std::endl;
00389 ofile << "POINTS " << n_nodes << " float" << std::endl;
00390 for (size_t i=0; i<n_nodes; ++i) { ofile << _8s()<<N[i].X() << _8s() << N[i].Y() << _8s() << N[i].Z() << std::endl; }
00391 ofile << std::endl;
00392 ofile << "CELLS "<< n_elems << " " << n_data << std::endl;
00393 for (size_t i=0; i<n_elems; ++i) { ofile << active_E[i]->nNodes() << " ";
00394 for (int j=0; j<active_E[i]->nNodes(); ++j) { ofile << active_E[i]->Connects()[j]->Number() << " "; }
00395 ofile << std::endl; }
00396 ofile << std::endl;
00397 ofile << "CELL_TYPES " << n_elems << std::endl;
00398 for (size_t i=0; i<n_elems; ++i) { ofile << active_E[i]->VTKCellType() << std::endl; }
00399 ofile << std::endl;
00400 ofile << "POINT_DATA " << n_nodes << std::endl;
00401
00402 { ofile << "VECTORS " << "Disp float" << std::endl;
00403 for (size_t j=0; j<n_nodes; ++j)
00404 if (N[j].HasVar(DUX) && N[j].HasVar(DUY) && N[j].HasVar(DUZ))
00405 { ofile << _8s() << N[j].DOFVar(DUX).EssentialVal
00406 << _8s() << N[j].DOFVar(DUY).EssentialVal
00407 << _8s() << N[j].DOFVar(DUZ).EssentialVal << std::endl; }
00408 else { ofile << _8s() << 0 << _8s() << 0 << _8s() << 0 << std::endl; }
00409 ofile << std::endl; }
00410
00411 for (size_t i=0; i<n_comps; i++) { ofile << "SCALARS " << Labels[i] << " float 1" << std::endl;
00412 ofile << "LOOKUP_TABLE default" << std::endl;
00413 for (size_t j=0; j<n_nodes; ++j) { ofile << _8s() << DataValues(j,i) << std::endl; }
00414 ofile << std::endl; }
00415
00416 ofile.close();
00417
00418 }
00419
00420 inline void Output::OpenFiles()
00421 {
00422
00423 for (int i_node=0; i_node<_n_out_nodes; ++i_node)
00424 {
00425
00426 _a_node_files[i_node].open(_idata->fnOutNODES[i_node].GetSTL().c_str(), std::ios::out);
00427 if (_a_node_files[i_node].fail())
00428 throw new Fatal(_("Output::Initialize: Could not open file < %s >"), _idata->fnOutNODES[i_node].c_str());
00429 }
00430
00431
00432 for (int i_elem=0; i_elem<_n_out_elems; ++i_elem)
00433 {
00434
00435 _a_elem_files[i_elem].open(_idata->fnOutELEMS[i_elem].GetSTL().c_str(), std::ios::out);
00436 if (_a_elem_files[i_elem].fail())
00437 throw new Fatal(_("Output::Initialize: Could not open file < %s >"), _idata->fnOutELEMS[i_elem].c_str());
00438 }
00439 }
00440
00441 inline void Output::CloseFiles()
00442 {
00443
00444 for (int i_node=0; i_node<_n_out_nodes; ++i_node)
00445 _a_node_files[i_node].close();
00446
00447
00448 for (int i_elem=0; i_elem<_n_out_elems; ++i_elem)
00449 _a_elem_files[i_elem].close();
00450 }
00451
00452 inline void Output::VTKCreateFile(String const & Filename)
00453 {
00454
00455 std::ofstream of(Filename.GetSTL().c_str(), std::ios::out);
00456
00457
00458 Array<FEM::Node> const & N = _data->Nodes;
00459 Array<FEM::Element*> const & E = _data->Elements;
00460
00461
00462 int n_nodes = N.size();
00463 int n_elems = E.size();
00464
00465
00466 int n_data = 0;
00467 for (int i=0; i<n_elems; ++i) n_data += 1 + E[i]->nNodes();
00468
00469
00470 of << "# vtk DataFile Version 3.0" << "\n";
00471 of << "MechSys/FEM - "<<_idata->projNAME << "\n";
00472 of << "ASCII" << "\n";
00473 of << "DATASET UNSTRUCTURED_GRID" << "\n";
00474 of << "\n";
00475 of << "POINTS "<<n_nodes<<" float" << "\n"; for (int i=0; i<n_nodes; ++i) {
00476 of << _8s()<<N[i].X() <<_8s()<<N[i].Y() <<_8s()<<N[i].Z() << "\n"; }
00477 of << "\n";
00478 of << "CELLS "<<n_elems<<" "<<n_data << "\n"; for (int i=0; i<n_elems; ++i) {
00479 of << E[i]->nNodes()<< " "; for (int j=0; j<E[i]->nNodes(); ++j)
00480 { of<<E[i]->Connects()[j]->Number()<< " "; } of << "\n"; }
00481 of << "\n";
00482 of << "CELL_TYPES " << n_elems << "\n"; for (int i=0; i<n_elems; ++i) {
00483 of << E[i]->VTKCellType() << "\n"; }
00484 of << "\n";
00485 of << "POINT_DATA " << n_nodes << "\n";
00486 of << "SCALARS Zcoordinate float 1" << "\n";
00487 of << "LOOKUP_TABLE default" << "\n"; for (int i=0; i<n_nodes; ++i) {
00488 of << _8s()<<N[i].Z() << "\n"; }
00489
00490
00491 of.close();
00492
00493 }
00494
00495 inline void Output::VTKAppendData(String const & Filename, int iStage, int iInc)
00496 {
00497
00498 std::ofstream of(Filename.GetSTL().c_str(), std::ios::app);
00499
00500
00501
00502 Array<FEM::Element*> const & E = _data->Elements;
00503
00504
00505
00506 int n_elems = E.size();
00507
00508
00509 of << "\n";
00510 of << "CELL_DATA " << n_elems << "\n"; if (_idata->vtkOutScalar1) {
00511 of << "SCALARS "<<iStage<<"-"<<iInc<<"-Scalar1 float 1" << "\n";
00512 of << "LOOKUP_TABLE default" << "\n"; for (int i=0; i<n_elems; ++i) {
00513 of << _8s()<<E[i]->OutScalar1() << "\n"; } } if (_idata->vtkOutScalar2) {
00514 of << "SCALARS "<<iStage<<"-"<<iInc<<"-Scalar2 float 1" << "\n";
00515 of << "LOOKUP_TABLE default" << "\n"; for (int i=0; i<n_elems; ++i) {
00516 of << _8s()<<E[i]->OutScalar2() << "\n"; } } if (_idata->vtkOutTensor1) {
00517 of << "TENSORS "<<iStage<<"-"<<iInc<<"-Tensor1 float" << "\n"; for (int i=0; i<n_elems; ++i) {
00518 String ten; E[i]->OutTensor1(ten); of << ten << "\n"; } } if (_idata->vtkOutTensor2) {
00519 of << "TENSORS "<<iStage<<"-"<<iInc<<"-Tensor2 float" << "\n"; for (int i=0; i<n_elems; ++i) {
00520 String ten; E[i]->OutTensor2(ten); of << ten << "\n"; } } if (_idata->vtkOutTensor3) {
00521 of << "TENSORS "<<iStage<<"-"<<iInc<<"-Tensor3 float" << "\n"; for (int i=0; i<n_elems; ++i) {
00522 String ten; E[i]->OutTensor3(ten); of << ten << "\n"; } }
00523 of << "\n";
00524
00525
00526
00527
00528
00529
00530
00531 of.close();
00532
00533 }
00534
00535 };
00536
00537 #endif // MECHSYS_FEM_OUTPUT_H
00538
00539