1 /** @file gsWriteParaviewMultiPhysics.cpp
2
3 @brief Provides implementation for gsWriteParaviewMultiPhysics.h
4
5 This file is part of the G+Smo library.
6
7 This Source Code Form is subject to the terms of the Mozilla Public
8 License, v. 2.0. If a copy of the MPL was not distributed with this
9 file, You can obtain one at http://mozilla.org/MPL/2.0/.
10
11 Author(s): A. Shamanskiy (TU Kaiserslautern)
12 Inspired by gsWriteParaview.hpp by A. Mantzaflaris
13 */
14
15 #include <gsElasticity/gsWriteParaviewMultiPhysics.h>
16 #include <gsUtils/gsPointGrid.h>
17 #include <gsUtils/gsMesh/gsMesh.h>
18 #include <gsCore/gsFunction.h>
19 #include <gsCore/gsField.h>
20 #include <gsIO/gsWriteParaview.h>
21 #include <gsElasticity/gsGeoUtils.h>
22
23
24 #define PLOT_PRECISION 11
25
26
27 namespace gismo
28 {
29
30
31 //---------- START REPEATED from gsWriteParaview.hpp
32
33 template<class T>
writeSingleControlNet(const gsGeometry<T> & Geo,std::string const & fn)34 void writeSingleControlNet(const gsGeometry<T> & Geo,
35 std::string const & fn)
36 {
37 const short_t d = Geo.parDim();
38 gsMesh<T> msh;
39 Geo.controlNet(msh);
40 const short_t n = Geo.geoDim();
41 if ( n == 1 )
42 {
43 gsMatrix<T> anch = Geo.basis().anchors();
44 // Lift vertices at anchor positions
45 for (std::size_t i = 0; i!= msh.numVertices(); ++i)
46 {
47 msh.vertex(i)[d] = msh.vertex(i)[0];
48 msh.vertex(i).topRows(d) = anch.col(i);
49 }
50 }
51 else if (n>3)
52 {
53 gsDebug<<"Writing 4th coordinate\n";
54 const gsMatrix<T> & cp = Geo.coefs();
55 gsWriteParaviewPoints<T>(cp.transpose(), fn );
56 return;
57 }
58
59 gsWriteParaview(msh, fn, false);
60 }
61
62 template<class T>
writeSingleCompMesh(const gsBasis<T> & basis,const gsGeometry<T> & Geo,std::string const & fn,unsigned resolution=8)63 void writeSingleCompMesh(const gsBasis<T> & basis, const gsGeometry<T> & Geo,
64 std::string const & fn, unsigned resolution = 8)
65 {
66 gsMesh<T> msh(basis, resolution);
67 Geo.evaluateMesh(msh);
68 gsWriteParaview(msh, fn, false);
69 }
70
71
72 //---------- END REPEATED from gsWriteParaview.hpp
73
74
75 template<class T>
gsWriteParaviewMultiPhysics(std::map<std::string,const gsField<T> * > fields,std::string const & fn,unsigned npts,bool mesh,bool ctrlNet)76 void gsWriteParaviewMultiPhysics(std::map<std::string, const gsField<T>*> fields,
77 std::string const & fn,
78 unsigned npts, bool mesh, bool ctrlNet)
79 {
80 const unsigned numP = fields.begin()->second->patches().nPatches();
81 gsParaviewCollection collection(fn);
82 std::string fileName = fn.substr(fn.find_last_of("/\\")+1); // file name without a path
83
84 for ( unsigned i=0; i < numP; ++i )
85 {
86 const gsBasis<> & dom = fields.begin()->second->isParametrized() ?
87 fields.begin()->second->igaFunction(i).basis() : fields.begin()->second->patch(i).basis();
88
89 gsWriteParaviewMultiPhysicsSinglePatch( fields, i, fn + util::to_string(i), npts);
90 collection.addPart(fileName + util::to_string(i), ".vts");
91
92 if ( mesh )
93 {
94 writeSingleCompMesh(dom, fields.begin()->second->patch(i), fn + util::to_string(i) + "_mesh");
95 collection.addPart(fileName + util::to_string(i) + "_mesh", ".vtp");
96 }
97 if ( ctrlNet ) // Output the control net
98 {
99 writeSingleControlNet(fields.begin()->second->patch(i), fn + util::to_string(i) + "_cnet");
100 collection.addPart(fileName + util::to_string(i) + "_cnet", ".vtp");
101 }
102
103 }
104 collection.save();
105 }
106
107 template<class T>
gsWriteParaviewMultiPhysicsTimeStep(std::map<std::string,const gsField<T> * > fields,std::string const & fn,gsParaviewCollection & collection,int time,unsigned npts)108 void gsWriteParaviewMultiPhysicsTimeStep(std::map<std::string, const gsField<T> *> fields, std::string const & fn,
109 gsParaviewCollection & collection, int time, unsigned npts)
110 {
111 const unsigned numP = fields.begin()->second->patches().nPatches();
112 std::string fileName = fn.substr(fn.find_last_of("/\\")+1); // file name without a path
113
114 for ( size_t p = 0; p < numP; ++p)
115 {
116 gsWriteParaviewMultiPhysicsSinglePatch(fields,p,fn + util::to_string(time) + "_" + util::to_string(p),npts);
117 collection.addTimestep(fileName + util::to_string(time) + "_",p,time,".vts");
118 }
119
120 }
121
122 template<class T>
gsWriteParaviewMultiPhysicsSinglePatch(std::map<std::string,const gsField<T> * > fields,const unsigned patchNum,std::string const & fn,unsigned npts)123 void gsWriteParaviewMultiPhysicsSinglePatch(std::map<std::string,const gsField<T> *> fields,
124 const unsigned patchNum,
125 std::string const & fn,
126 unsigned npts)
127 {
128 const gsGeometry<> & geometry = fields.begin()->second->patches().patch(patchNum);
129 const short_t n = geometry.targetDim();
130 const short_t d = geometry.domainDim();
131
132 gsMatrix<> ab = geometry.support();
133 gsVector<> a = ab.col(0);
134 gsVector<> b = ab.col(1);
135 gsVector<unsigned> np = distributePoints<T>(geometry,npts);
136 gsMatrix<> pts = gsPointGrid(a,b,np);
137
138 gsMatrix<> eval_geo = geometry.eval(pts);
139 std::map<std::string, gsMatrix<> > data;
140 for (typename std::map<std::string,const gsField<T> *>::iterator it = fields.begin(); it != fields.end(); it++)
141 {
142 data[it->first] = it->second->isParametric() ?
143 it->second->function(patchNum).eval(pts) : it->second->function(patchNum).eval(eval_geo);
144
145 if ( data[it->first].rows() == 2 )
146 {
147 data[it->first].conservativeResize(3,eval_geo.cols() );
148 data[it->first].row(2).setZero();
149 }
150 }
151
152 if (3 -d > 0)
153 {
154 np.conservativeResize(3);
155 np.bottomRows(3-d).setOnes();
156 }
157 else if (d > 3)
158 {
159 gsWarn<< "Cannot plot 4D data.\n";
160 return;
161 }
162
163 if ( 3 - n > 0 )
164 {
165 eval_geo.conservativeResize(3,eval_geo.cols() );
166 eval_geo.bottomRows(3-n).setZero();
167 }
168 else if (n > 3)
169 {
170 gsWarn<< "Data is more than 3 dimensions.\n";
171 }
172
173 /*for (typename std::map<std::string, gsMatrix<> >::iterator it = data.begin(); it != data.end(); it++)
174 {
175 if ( it->second.rows() > 1 )
176 {
177 it->second.conservativeResize(3,eval_geo.cols() );
178 it->second.bottomRows( 3-dd ).setZero();
179 }
180 }*/
181
182 gsWriteParaviewMultiTPgrid(eval_geo, data, np.template cast<index_t>(), fn);
183 }
184
185 template<class T>
gsWriteParaviewMultiTPgrid(gsMatrix<T> const & points,std::map<std::string,gsMatrix<T>> & data,const gsVector<index_t> & np,std::string const & fn)186 void gsWriteParaviewMultiTPgrid(gsMatrix<T> const& points,
187 std::map<std::string, gsMatrix<T> >& data,
188 const gsVector<index_t> & np,
189 std::string const & fn)
190 {
191 const int n = points.rows();
192
193 std::string mfn(fn);
194 mfn.append(".vts");
195 std::ofstream file(mfn.c_str());
196 file << std::fixed; // no exponents
197 file << std::setprecision (PLOT_PRECISION);
198
199 file <<"<?xml version=\"1.0\"?>\n";
200 file <<"<VTKFile type=\"StructuredGrid\" version=\"0.1\">\n";
201 file <<"<StructuredGrid WholeExtent=\"0 "<< np(0)-1<<" 0 "<<np(1)-1<<" 0 "
202 << (np.size()>2 ? np(2)-1 : 0) <<"\">\n";
203 file <<"<Piece Extent=\"0 "<< np(0)-1<<" 0 "<<np(1)-1<<" 0 "
204 << (np.size()>2 ? np(2)-1 : 0) <<"\">\n";
205
206 file <<"<PointData>\n";
207 for (typename std::map<std::string, gsMatrix<T> >::iterator it = data.begin(); it != data.end(); it++)
208 {
209 file <<"<DataArray type=\"Float32\" Name=\""<< it->first <<"\" format=\"ascii\" NumberOfComponents=\""<< ( it->second.rows()==1 ? 1 : 3) <<"\">\n";
210 if ( it->second.rows()==1 )
211 for ( index_t j=0; j<it->second.cols(); ++j)
212 file<< it->second.at(j) <<" ";
213 else
214 {
215 for ( index_t j=0; j<it->second.cols(); ++j)
216 {
217 for ( index_t i=0; i!=it->second.rows(); ++i)
218 file<< it->second(i,j) <<" ";
219 for ( index_t i=it->second.rows(); i<3; ++i)
220 file<<"0 ";
221 }
222 }
223 file <<"</DataArray>\n";
224 }
225 file <<"</PointData>\n";
226 file <<"<Points>\n";
227 file <<"<DataArray type=\"Float32\" NumberOfComponents=\"3\">\n";
228 for ( index_t j=0; j<points.cols(); ++j)
229 {
230 for ( index_t i=0; i!=n; ++i)
231 file<< points(i,j) <<" ";
232 for ( index_t i=n; i<3; ++i)
233 file<<"0 ";
234 }
235 file <<"</DataArray>\n";
236 file <<"</Points>\n";
237 file <<"</Piece>\n";
238 file <<"</StructuredGrid>\n";
239 file <<"</VTKFile>\n";
240
241 file.close();
242 }
243
244
245 }
246
247 #undef PLOT_PRECISION
248