00001 #include "gnuplotdeal.h"
00002 #include <numerics/data_out.h>
00003 #include <iostream>
00004 #include <fstream>
00005 #include <math.h>
00006 #include <lac/full_matrix.h>
00007 #include <base/quadrature_lib.h>
00008 #include <grid/grid_out.h>
00009 using std::endl;
00010 using std::cout;
00011
00012
00021 void GnuPlotDeal::plotDealSolution(DoFHandler<1> &dof,const VecDouble &sol,std::string name,std::string strComp)
00022 {
00023 string strFileName = getDataFileName();
00024 m_List.push_back(strFileName);
00025 std::ofstream fout(strFileName.c_str());
00026 if (m_bPlot3D)
00027 {
00028 std::cerr << "Graficos 2D e 3D na mesma cena nao pode ser feito\n";
00029 this->finalize();
00030 }
00031 DataOut<1> data_out;
00032 data_out.attach_dof_handler (dof);
00033 data_out.add_data_vector ((DealVecDouble&)sol, name);
00034 data_out.build_patches ();
00035 data_out.write_gnuplot (fout);
00036 if (!m_bPlot) {
00037 fprintf(m_file,"plot \"%s\" %s title \"%s\" ",strFileName.c_str(),strComp.c_str(),name.c_str());
00038 m_bPlot = true;
00039 }
00040 else {
00041 fprintf(m_file,", \"%s\" %s title \"%s\"",strFileName.c_str(),strComp.c_str(),name.c_str());
00042 }
00043 return;
00044
00045
00046 }
00047
00048
00049
00058 void GnuPlotDeal::plotDealSolutionDegAtVertex(DoFHandler<1> &dofs,VecDouble &sol,vector<int> degOfVertex,std::string name,std::string strComp)
00059 {
00060 string strFileName = getDataFileName();
00061 m_List.push_back(strFileName);
00062 std::ofstream fout(strFileName.c_str());
00063 if (m_bPlot3D)
00064 {
00065 std::cerr << "Graficos 2D e 3D na mesma cena nao pode ser feito\n";
00066 this->finalize();
00067 }
00068 fout << "#Plotted For Discontinuos Galerkin" << endl;
00069 const unsigned int dofs_per_cell = dofs.get_fe().dofs_per_cell;
00070 std::vector<unsigned int> local_dof_indices (dofs_per_cell);
00071
00072 DoFHandler<1>::active_cell_iterator cell = dofs.begin_active(),
00073 endc = dofs.end();
00074 for (; cell!=endc; ++cell)
00075 {
00076 cell->get_dof_indices (local_dof_indices);
00077 for (unsigned i=0;i<2;i++)
00078 fout << cell->vertex(i)[0] << " " << sol(cell->vertex_dof_index(i,degOfVertex[i])) << endl;
00079 fout << endl;
00080 }
00081 fout.close();
00082 if (!m_bPlot) {
00083 fprintf(m_file,"plot \"%s\" %s title \"%s\" ",strFileName.c_str(),strComp.c_str(),name.c_str());
00084 m_bPlot = true;
00085 }
00086 else {
00087 fprintf(m_file,", \"%s\" %s title \"%s\"",strFileName.c_str(),strComp.c_str(),name.c_str());
00088 }
00089 return;
00090
00091
00092 }
00093
00094
00095
00096
00097
00104 void GnuPlotDeal::plotDealSolution(DoFHandler<2> &dofs,const VecDouble &sol,std::string name,std::string strComp)
00105 {
00106 unsigned i=0,deg=0;
00107 string strFileName = getDataFileName();
00108 m_List.push_back(strFileName);
00109 std::ofstream fout(strFileName.c_str());
00110 fout << "#Plotando x y + (todos graus de liberdade) " << deg << std::endl;
00111 if (m_bPlot)
00112 {
00113 std::cerr << "Graficos 2D e 3D na mesma cena nao pode ser feito\n";
00114 this->finalize();
00115 }
00116 const unsigned int dofs_per_cell = dofs.get_fe().dofs_per_cell;
00117 std::vector<unsigned int> local_dof_indices (dofs_per_cell);
00118 DoFHandler<2>::active_cell_iterator cell = dofs.begin_active(),
00119 endc = dofs.end();
00120 unsigned n_components = dofs.get_fe().n_components();
00121 for (; cell!=endc; ++cell)
00122 {
00123 cell->get_dof_indices (local_dof_indices);
00124 for (i=0;i<GeometryInfo<2>::vertices_per_cell;i++)
00125 {
00126 fout << cell->vertex(i)[0] << " " << cell->vertex(i)[1] << " ";
00127 for (deg=0;deg<n_components;deg++)
00128 fout << sol(cell->vertex_dof_index(i,deg)) << " ";
00129 fout << endl;
00130 }
00131 i=0;
00132 fout << cell->vertex(i)[0] << " " << cell->vertex(i)[1] << " ";
00133 for (deg=0;deg<n_components;deg++)
00134 fout << sol(cell->vertex_dof_index(i,deg)) << " ";
00135 fout << endl<<endl<<endl;
00136 }
00137
00138 if (!m_bPlot3D) {
00139 fprintf(m_file,"splot \"%s\" %s title \"%s\" ",strFileName.c_str(),strComp.c_str(),name.c_str());
00140 m_bPlot3D = true;
00141 }
00142 else {
00143 fprintf(m_file,", \"%s\" %s title \"%s\"",strFileName.c_str(),strComp.c_str(),name.c_str());
00144 }
00145
00146 return;
00147 }
00148
00149
00160 void GnuPlotDeal::plotDealSolution(DoFHandler<2> &dofs,const VecDouble &sol,std::string name,std::string strComp,unsigned deg)
00161 {
00162 unsigned i=0;
00163 string strFileName = getDataFileName();
00164 m_List.push_back(strFileName);
00165 std::ofstream fout(strFileName.c_str());
00166 fout << "#Plotando x y f(x,y) no grau de liberdade " << deg << std::endl;
00167
00168 const unsigned int dofs_per_cell = dofs.get_fe().dofs_per_cell;
00169 std::vector<unsigned int> local_dof_indices (dofs_per_cell);
00170
00171 DoFHandler<2>::active_cell_iterator cell = dofs.begin_active(),
00172 endc = dofs.end();
00173 fout.precision(10);
00174 for (; cell!=endc; ++cell)
00175 {
00176 cell->get_dof_indices (local_dof_indices);
00177 for (i=0;i<GeometryInfo<2>::vertices_per_cell;i++)
00178 fout << cell->vertex(i)[0] << " " << cell->vertex(i)[1] << " " << sol(cell->vertex_dof_index(i,deg)) << endl;
00179 i=0;
00180 fout << cell->vertex(i)[0] << " " << cell->vertex(i)[1] << " " << sol(cell->vertex_dof_index(i,deg)) << endl << endl << endl;
00181 }
00182 fout.close();
00183
00184 plotPoints3d(strFileName,name,strComp);
00185
00186 return;
00187 }
00188
00189
00190
00191
00192
00207 void GnuPlotDeal::plotDealVectorFieldAtVertices(DoFHandler<2> &dofs,VecDouble &sol,unsigned deg1,unsigned deg2,bool bScale,double factor,bool bNormalize, std::string name,std::string strComp)
00208 {
00209 std::vector<Point<2> > Verts(4);
00210 Verts.push_back(Point<2> (0,0));
00211 Verts.push_back(Point<2> (1,0));
00212 Verts.push_back(Point<2> (1,1));
00213 Verts.push_back(Point<2> (0,1));
00214
00215 Quadrature<2> quadrature(Verts);
00216 FEValues<2> fe_values (dofs.get_fe(), quadrature,
00217 UpdateFlags(update_values|update_q_points));
00218 std::ofstream fout;
00219 std::string strFileName = newDataFile(fout);
00220 unsigned dofs_per_cell = dofs.get_fe().dofs_per_cell;
00221 unsigned n_q_points = fe_values.n_quadrature_points;
00222 vector<unsigned> dofIndices(dofs.get_fe().dofs_per_cell);
00223
00224 DoFHandler<2>::active_cell_iterator cell = dofs.begin_active(),
00225 endc = dofs.end();
00226 Point<2> Q,V,QAux,VAux;
00227 double maxHyp=0.0;
00228 for (;cell!=endc; ++cell)
00229 {
00230 fe_values.reinit(cell);
00231 cell->get_dof_indices(dofIndices);
00232 for (unsigned q=0;q<n_q_points;q++)
00233 {
00234 Q = fe_values.quadrature_point(q);
00235 V*=0;
00236 for (unsigned i=0;i<dofs_per_cell;i++)
00237 {
00238 V(0) += sol(dofIndices[i])*fe_values.shape_value_component(i,q,deg1);
00239 V(1) += sol(dofIndices[i])*fe_values.shape_value_component(i,q,deg2);
00240 }
00241 double hyp = hypot(V(0),V(1));
00242 if (bNormalize)
00243 {
00244 V/=hyp;
00245 hyp=1.0;
00246 }
00247 if (hyp > maxHyp) maxHyp = hyp;
00248
00249 fout << Q(0) << " " << Q(1) << " " << V(0) << " " << V(1) << " " << endl;
00250 if (q == 0)
00251 {
00252 QAux = Q;
00253 VAux = V;
00254 }
00255 }
00256 fout << QAux(0) << " " << QAux(1) << " " << VAux(0) << " " << VAux(1) << " " << endl;
00257 fout << endl <<endl;
00258 }
00259 fout.close();
00260 double rad = dofs.begin_active()->diameter()/2.0;
00261
00262 char strAux[100];
00263 if (bScale)
00264 sprintf(strAux," using 1:2:(%g*$3):(%g*$4) with vector ",factor*rad/maxHyp,factor*rad/maxHyp);
00265 else
00266 sprintf(strAux," with vector ");
00267
00268 std::string strWith = strAux;
00269 strWith+=strComp;
00270 plotPoints(strFileName,name,strWith);
00271 }
00272
00273
00274
00275
00276
00277
00278
00279
00286 void GnuPlotDeal::plotDealVectorField2D(DoFHandler<2> &dofs,VecDouble &sol1,VecDouble &sol2,std::string name,bool bScale,std::string strComp)
00287 {
00288
00289
00290 string strFileName = getDataFileName();
00291 double aux,MaxNorm = 0.0,MaxDiameter=0.0;;
00292 int NumPoints = 0;
00293 m_List.push_back(strFileName);
00294 std::ofstream fout(strFileName.c_str());
00295 if (m_bPlot3D)
00296 {
00297 std::cerr << "Graficos 2D e 3D na mesma cena nao pode ser feito\n";
00298 this->finalize();
00299 }
00300 fout << "#Vector Field Plotted by Deal. Format <X,Y,DX,DY>\n";
00301 const unsigned int dofs_per_cell = dofs.get_fe().dofs_per_cell;
00302
00303 std::vector<unsigned int> local_dof_indices (dofs_per_cell);
00304
00305 DoFHandler<2>::active_cell_iterator cell = dofs.begin_active(),
00306 endc = dofs.end();
00307 for (; cell!=endc; ++cell)
00308 {
00309 cell->get_dof_indices (local_dof_indices);
00310 Point2D p;
00311 double diam = cell->diameter();
00312 if (diam > MaxDiameter)
00313 {
00314 MaxDiameter = diam;
00315 }
00316 for (unsigned int i=0; i<dofs_per_cell; ++i,NumPoints++)
00317 {
00318 p(0) = sol1(local_dof_indices[i]);
00319 p(1) = sol2(local_dof_indices[i]);
00320 aux = sqrt(p.square());
00321 if (aux > MaxNorm)
00322 MaxNorm = aux;
00323 }
00324 }
00325 FullMatrix<double> M(NumPoints,4);
00326 int lineCnt = 0;
00327 double factor = MaxDiameter/(MaxNorm*2.0);
00328
00329 if (bScale)
00330 {
00331 cell = dofs.begin_active();
00332 for (; cell!=endc; ++cell)
00333 {
00334 cell->get_dof_indices (local_dof_indices);
00335 for (unsigned int i=0; i<dofs_per_cell; ++i,lineCnt++)
00336 {
00337 M(lineCnt,0) = cell->vertex(i)[0];
00338 M(lineCnt,1) = cell->vertex(i)[1];
00339 M(lineCnt,2) = sol1(local_dof_indices[i])*factor;
00340 M(lineCnt,3) = sol2(local_dof_indices[i])*factor;
00341 }
00342 }
00343 }
00344 else{
00345 cell = dofs.begin_active();
00346 for (; cell!=endc; ++cell)
00347 {
00348 cell->get_dof_indices (local_dof_indices);
00349 for (unsigned int i=0; i<dofs_per_cell; ++i,lineCnt++)
00350 {
00351 M(lineCnt,0) = cell->vertex(i)[0];
00352 M(lineCnt,1) = cell->vertex(i)[1];
00353 M(lineCnt,2) = sol1(local_dof_indices[i]);
00354 M(lineCnt,3) = sol2(local_dof_indices[i]);
00355 factor = 1.0/hypot(M(lineCnt,2),M(lineCnt,3));
00356 M(lineCnt,2)*=factor*MaxDiameter/2.0;
00357 M(lineCnt,3)*=factor*MaxDiameter/2.0;
00358 }
00359 }
00360 }
00361
00362 M.print_formatted(fout,6,false,0,"0,0");
00363 if (!m_bPlot) {
00364 fprintf(m_file,"plot \"%s\" with vector %s title \"%s\" ",strFileName.c_str(),strComp.c_str(),name.c_str());
00365 m_bPlot = true;
00366 }
00367 else {
00368 fprintf(m_file,", \"%s\" with vector %s title \"%s\"",strFileName.c_str(),strComp.c_str(),name.c_str());
00369 }
00370 return;
00371 }
00372
00373
00386 void GnuPlotDeal::plotDealGradientField2D(DoFHandler<2> &dofs,VecDouble &sol,unsigned deg,bool bScale,double factor,bool bNormalize,std::string name,std::string strComp)
00387 {
00388 std::vector<Point<2> > Verts(4);
00389 Verts.push_back(Point<2> (0,0));
00390 Verts.push_back(Point<2> (1,0));
00391 Verts.push_back(Point<2> (1,1));
00392 Verts.push_back(Point<2> (0,1));
00393
00394 Quadrature<2> quadrature(Verts);
00395 FEValues<2> fe_values (dofs.get_fe(), quadrature,
00396 UpdateFlags(update_gradients|update_q_points));
00397 std::ofstream fout;
00398 std::string strFileName = newDataFile(fout);
00399 unsigned dofs_per_cell = dofs.get_fe().dofs_per_cell;
00400 unsigned n_q_points = fe_values.n_quadrature_points;
00401 vector<unsigned> dofIndices(dofs.get_fe().dofs_per_cell);
00402
00403 DoFHandler<2>::active_cell_iterator cell = dofs.begin_active(),
00404 endc = dofs.end();
00405 Point<2> Q,QAux;
00406 Tensor<1,2> V,VAux;
00407 double maxHyp=0.0;
00408 for (;cell!=endc; ++cell)
00409 {
00410 fe_values.reinit(cell);
00411 cell->get_dof_indices(dofIndices);
00412 for (unsigned q=0;q<n_q_points;q++)
00413 {
00414 Q = fe_values.quadrature_point(q);
00415 V*=0;
00416 for (unsigned i=0;i<dofs_per_cell;i++)
00417 {
00418 V += sol(dofIndices[i])*fe_values.shape_grad_component(i,q,deg);
00419 }
00420 double hyp = hypot(V[0],V[1]);
00421 if (bNormalize)
00422 {
00423 V/=hyp;
00424 hyp=1.0;
00425 }
00426 if (hyp > maxHyp) maxHyp = hyp;
00427
00428 fout << Q(0) << " " << Q(1) << " " << V[0] << " " << V[1] << " " << endl;
00429 if (q == 0)
00430 {
00431 QAux = Q;
00432 VAux = V;
00433 }
00434 }
00435 fout << QAux(0) << " " << QAux(1) << " " << VAux[0] << " " << VAux[1] << " " << endl;
00436 fout << endl <<endl;
00437 }
00438 fout.close();
00439 double rad = dofs.begin_active()->diameter()/2.0;
00440
00441 char strAux[100];
00442 if (bScale)
00443 sprintf(strAux," using 1:2:(%g*$3):(%g*$4) with vector ",factor*rad/maxHyp,factor*rad/maxHyp);
00444 else
00445 sprintf(strAux," with vector ");
00446
00447 std::string strWith = strAux;
00448 strWith+=strComp;
00449 plotPoints(strFileName,name,strWith);
00450 }
00451
00452
00453
00454
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00475 void GnuPlotDeal::plotDealGrid2D(const Triangulation<2> &tri,std::string name,std::string strComp)
00476 {
00477 string strFileName = getDataFileName();
00478 m_List.push_back(strFileName);
00479 std::ofstream fout(strFileName.c_str());
00480 GridOut grid_out;
00481 grid_out.write_gnuplot(tri, fout);
00482 this->plotPoints(strFileName,name,strComp);
00483
00484 }
00485
00486
00493 void GnuPlotDeal::plotDealVectorField2DGaussPoints(DoFHandler<2> &dofs,VecDouble &sol1,VecDouble &sol2,std::string name,bool bScale,double factor,std::string strComp)
00494 {
00495 string strFileName = getDataFileName();
00496 m_List.push_back(strFileName);
00497 std::ofstream fout(strFileName.c_str());
00498 double norm,MaxNorm=0.0,MaxMeasure=0.0;
00499 unsigned k;
00500 if (m_bPlot3D)
00501 {
00502 std::cerr << "Graficos 2D e 3D na mesma cena nao pode ser feito\n";
00503 this->finalize();
00504 }
00505
00506 fout << "#Plot Gradients on Gauss Points";
00507 QGauss<2> quadrature_formula(2);
00508 FEValues<2> fe_values (dofs.get_fe(), quadrature_formula,
00509 UpdateFlags(update_values |
00510 update_gradients | update_q_points |
00511 update_JxW_values));
00512
00513 const unsigned int dofs_per_cell = dofs.get_fe().dofs_per_cell;
00514 const unsigned int n_q_points = quadrature_formula.n_quadrature_points;
00515
00516 FullMatrix<double> M(dofs.get_tria().n_active_quads()*4,4);
00517 M = 0;
00518 std::vector<unsigned int> local_dof_indices (dofs_per_cell);
00519
00520 DoFHandler<2>::active_cell_iterator cell = dofs.begin_active(),
00521 endc = dofs.end();
00522 for (k=0; cell!=endc; ++cell)
00523 {
00524 fe_values.reinit (cell);
00525 cell->get_dof_indices(local_dof_indices);
00526 for (unsigned int q_point=0; q_point<n_q_points; ++q_point,k++)
00527 {
00528 M(k,0) = (fe_values.quadrature_point(q_point)[0]);
00529 M(k,1) = (fe_values.quadrature_point(q_point)[1]);
00530 for (unsigned int i=0; i<dofs_per_cell; ++i)
00531 {
00532 M(k,2)+= sol1(local_dof_indices[i])*fe_values.shape_value(i, q_point);
00533 M(k,3)+= sol2(local_dof_indices[i])*fe_values.shape_value(i, q_point);
00534 }
00535 double hyp = hypot(M(k,2),M(k,3));
00536 if (!bScale)
00537 {
00538 if (hyp >1.0E-08)
00539 {
00540 norm = 1.0/hyp * sqrt(fe_values.JxW (q_point))*factor;
00541 }
00542 else{
00543 norm = 0.0;
00544 }
00545 M(k,2)*=norm;
00546 M(k,3)*=norm;
00547 }
00548 else
00549 {
00550 if (hyp > MaxNorm) MaxNorm=hyp;
00551 if (MaxMeasure < fe_values.JxW(q_point)) MaxMeasure = fe_values.JxW(q_point);
00552 }
00553 }
00554 }
00555 if (bScale)
00556 {
00557 double c = sqrt(MaxMeasure)/(MaxNorm)*factor;
00558 for (k=0;k<M.size(0);k++)
00559 {
00560 M(k,2)*=c;
00561 M(k,3)*=c;
00562 }
00563 }
00564 M.print_formatted(fout,6,false,0,"0,0");
00565
00566 if (!m_bPlot) {
00567 fprintf(m_file,"plot \"%s\" with vector %s title \"%s\" ",strFileName.c_str(),strComp.c_str(),name.c_str());
00568 m_bPlot = true;
00569 }
00570 else {
00571 fprintf(m_file,", \"%s\" with vector %s title \"%s\"",strFileName.c_str(),strComp.c_str(),name.c_str());
00572 }
00573
00574
00575
00576 }
00577
00578
00579
00580
00590 void GnuPlotDeal::plotDealDeformationGrid2D(DoFHandler<2> &dofs,const VecDouble &sol,std::string name,std::string strComp)
00591 {
00592 string strFileName = getDataFileName();
00593 m_List.push_back(strFileName);
00594 std::ofstream fout(strFileName.c_str());
00595 if (m_bPlot3D)
00596 {
00597 std::cerr << "Graficos 2D e 3D na mesma cena nao pode ser feito\n";
00598 this->finalize();
00599 }
00600 fout << "#Grid Plotted by Deal. Format <X,Y>\n";
00601 const unsigned int dofs_per_cell = dofs.get_fe().dofs_per_cell;
00602 std::vector<unsigned int> local_dof_indices (dofs_per_cell);
00603
00604 DoFHandler<2>::active_cell_iterator cell = dofs.begin_active(),
00605 endc = dofs.end();
00606 unsigned iVertex;
00607
00608
00609
00610 for (; cell!=endc; ++cell)
00611 {
00612
00613 for (iVertex=0;iVertex<GeometryInfo<2>::vertices_per_cell;iVertex++)
00614 {
00615 for (unsigned deg=0;deg<2;deg++)
00616 {
00617 fout << cell->vertex(iVertex)[deg]+sol(cell->vertex_dof_index(iVertex,deg)) << " ";
00618 }
00619
00620 fout << endl;
00621 }
00622 iVertex=0;
00623 for (unsigned deg=0;deg<2;deg++)
00624 {
00625 fout << cell->vertex(iVertex)[deg]+sol(cell->vertex_dof_index(iVertex,deg)) << " ";
00626 }
00627 fout << endl;
00628 fout<< endl;
00629 }
00630 fout.close();
00631
00632 if (!m_bPlot) {
00633 fprintf(m_file,"plot \"%s\" %s title \"%s\" ",strFileName.c_str(),strComp.c_str(),name.c_str());
00634 m_bPlot = true;
00635 }
00636 else {
00637 fprintf(m_file,", \"%s\" %s title \"%s\"",strFileName.c_str(),strComp.c_str(),name.c_str());
00638 }
00639
00640 return;
00641 }
00642
00643
00644
00645
00646
00653 void GnuPlotDeal::plotDealBoundary(DoFHandler<2> &dofs,unsigned boundary_indicator,std::string name,std::string strComp)
00654 {
00655 const FiniteElement<2> &fe = dofs.get_fe();
00656 std::vector<unsigned int> face_dofs (fe.dofs_per_face,DoFHandler<2>::invalid_dof_index);
00657 std::vector<Point<2> > dof_locations (face_dofs.size(), Point<2>());
00658 std::vector<double> dof_values_scalar (fe.dofs_per_face);
00659 std::vector<Point<2-1> > unit_support_points = fe.get_unit_face_support_points();
00660 assert(unit_support_points.size() !=0);
00661 assert(fe.is_primitive());
00662 Quadrature<2-1> aux_quad (unit_support_points);
00663 FEFaceValues<2> fe_values (fe, aux_quad, update_q_points);
00664 string strFileName = getDataFileName();
00665 m_List.push_back(strFileName);
00666 std::ofstream fout(strFileName.c_str());
00667 if (m_bPlot3D)
00668 {
00669 std::cerr << "Graficos 2D e 3D na mesma cena nao pode ser feito\n";
00670 this->finalize();
00671 }
00672 fout << "#Boundary Plotted by Deal for GnuPlot. Format <X,Y> <X,Y>\n";
00673
00674 DoFHandler<2>::active_cell_iterator cell = dofs.begin_active(),endc=dofs.end();
00675 for (; cell!=endc; ++cell)
00676 {
00677 for (unsigned int face_no = 0; face_no < GeometryInfo<2>::faces_per_cell;++face_no)
00678 {
00679 DoFHandler<2>::face_iterator face = cell->face(face_no);
00680 if (face->boundary_indicator() == boundary_indicator)
00681 {
00682 for (unsigned i=0;i<2;i++)
00683 {
00684 fout << face->vertex(i)(0) << " " << face->vertex(i)(1) << std::endl;
00685 }
00686 }
00687 fout << endl;
00688 }
00689 }
00690 fout.close();
00691 if (!m_bPlot) {
00692 fprintf(m_file,"plot \"%s\" %s title \"%s\" ",strFileName.c_str(),strComp.c_str(),name.c_str());
00693 m_bPlot = true;
00694 }
00695 else {
00696 fprintf(m_file,", \"%s\" %s title \"%s\"",strFileName.c_str(),strComp.c_str(),name.c_str());
00697 }
00698 }
00699
00700
00707 void GnuPlotDeal::plotDealPoints(vector<Point<2> > &points,string name,string strComp)
00708 {
00709 string strFileName = getDataFileName();
00710 m_List.push_back(strFileName);
00711 std::ofstream fout(strFileName.c_str());
00712
00713 if (m_bPlot3D)
00714 {
00715 std::cerr << "Graficos 2D e 3D na mesma cena nao pode ser feito\n";
00716 this->finalize();
00717 }
00718 fout << "#Plot Points sequence by Deal for GnuPlot. Format <X,Y> <X,Y>\n";
00719 for (unsigned i=0;i<points.size();i++)
00720 {
00721 fout << points[i](0) << " " << points[i](1) << endl;
00722
00723 }
00724 fout.close();
00725 if (!m_bPlot) {
00726 fprintf(m_file,"plot \"%s\" %s title \"%s\" ",strFileName.c_str(),strComp.c_str(),name.c_str());
00727 m_bPlot = true;
00728 }
00729 else {
00730 fprintf(m_file,", \"%s\" %s title \"%s\"",strFileName.c_str(),strComp.c_str(),name.c_str());
00731 }
00732 }
00733
00734
00741 void GnuPlotDeal::plotDealRadDeformationPoints(DoFHandler<2> &dofs,VecDouble &sol,double degree,double _TOLERANCE,std::string name,std::string strComp)
00742 {
00743 string strFileName = getDataFileName();
00744 m_List.push_back(strFileName);
00745 std::ofstream fout(strFileName.c_str());
00746 if (m_bPlot3D)
00747 {
00748 std::cerr << "Graficos 2D e 3D na mesma cena nao pode ser feito\n";
00749 this->finalize();
00750 }
00751 fout << "#Grid Plotted by Deal. Format <Inital Raio,End Raio>\n";
00752 const unsigned int dofs_per_cell = dofs.get_fe().dofs_per_cell;
00753 std::vector<unsigned int> local_dof_indices (dofs_per_cell);
00754
00755 DoFHandler<2>::active_cell_iterator cell = dofs.begin_active(),
00756 endc = dofs.end();
00757 unsigned iVertex;
00758
00759 degree=degree*M_PI/180.0;
00760 _TOLERANCE=_TOLERANCE*M_PI/180.0;
00761
00762 assert(dofs.get_fe().n_dofs_per_vertex()==2);
00763
00764 for (; cell!=endc; ++cell)
00765 {
00766
00767
00768
00769 for (iVertex=0;iVertex<GeometryInfo<2>::vertices_per_cell;iVertex++)
00770 {
00771 const Point<2> &P = cell->vertex(iVertex);
00772 if (fabs(atan(P[1]/P[0]) - degree) < _TOLERANCE)
00773 {
00774 fout << hypot(P[0],P[1]) << " " << hypot(sol(cell->vertex_dof_index(iVertex,0)),sol(cell->vertex_dof_index(iVertex,1))) << "\n";
00775 }
00776 }
00777 }
00778 fout.close();
00779
00780 if (!m_bPlot) {
00781 fprintf(m_file,"plot \"%s\" %s title \"%s\" ",strFileName.c_str(),strComp.c_str(),name.c_str());
00782 m_bPlot = true;
00783 }
00784 else {
00785 fprintf(m_file,", \"%s\" %s title \"%s\"",strFileName.c_str(),strComp.c_str(),name.c_str());
00786 }
00787
00788 return;
00789
00790 }
00791
00792
00793
00794
00806 void GnuPlotDeal::plotDealVectorField2DAtCentralPoints(Triangulation<2> &tria,VecDouble &V1,VecDouble &V2,double factor,bool bScale,std::string name,std::string strComp)
00807 {
00808 Triangulation<2>::active_cell_iterator cell = tria.begin_active(),
00809 endc = tria.end();
00810
00811 std::ofstream fout;
00812 std::string strFileName = newDataFile(fout);
00813 double MaxNorm=0,hyp;
00814 double rad = cell->diameter()/2.0;
00815 rad = sqrt(cell->measure())/2.0;
00816 for (;cell!=endc; ++cell)
00817 {
00818 Point2D XM = cell->barycenter();
00819 fout << XM(0) << " "
00820 << XM(1) << " "
00821 << V1(cell->index())*factor << " "
00822 << V2(cell->index())*factor << "\n\n";
00823 hyp=hypot(V1(cell->index()),V2(cell->index()));
00824 if (MaxNorm < hyp) MaxNorm=hyp;
00825
00826 }
00827 fout << std::endl;
00828 if (MaxNorm==0)
00829 MaxNorm= 1.0;
00830 char strAux[100];
00831 if (bScale)
00832 sprintf(strAux," using 1:2:(%g*$3):(%g*$4) with vector ",rad/MaxNorm,rad/MaxNorm);
00833 else
00834 sprintf(strAux," with vector ");
00835
00836 std::string strWith = strAux;
00837 strWith+=strComp;
00838 plotPoints(strFileName,name,strWith);
00839
00840 }
00841
00842
00843
00844
00855 void GnuPlotDeal::plotDealCentralValuesSolution(const Triangulation<2> &trig,const VecDouble &cSol,std::string name,std::string strComp)
00856 {
00857 assert(cSol.size() == trig.n_active_cells());
00858 string strFileName = getDataFileName();
00859 m_List.push_back(strFileName);
00860 std::ofstream fout(strFileName.c_str());
00861
00862 fout << "#Plotted For Const Central Values" << endl;
00863
00864 Triangulation<2>::active_cell_iterator cell = trig.begin_active(),
00865 endc = trig.end();
00866 for (; cell!=endc; ++cell)
00867 {
00868 fout << cell->vertex(0)[0] << " " << cell->vertex(0)[1] << " " << cSol(cell->index()) << "\n";
00869 fout << cell->vertex(1)[0] << " " << cell->vertex(1)[1] << " " << cSol(cell->index()) << "\n";
00870 fout << cell->vertex(2)[0] << " " << cell->vertex(2)[1] << " " << cSol(cell->index()) << "\n";
00871 fout << cell->vertex(3)[0] << " " << cell->vertex(3)[1] << " " << cSol(cell->index()) << "\n";
00872 fout << cell->vertex(0)[0] << " " << cell->vertex(0)[1] << " " << cSol(cell->index()) << "\n";
00873 fout << endl << endl;
00874 }
00875 fout.close();
00876
00877 plotPoints3d(strFileName,name,strComp);
00878
00879 return;
00880
00881 }
00882
00883 void GnuPlotDeal::plotDealSolutionAtVertices(Triangulation<2> &tria,const VecDouble &vSol,std::string name,std::string strComp)
00884 {
00885 assert(vSol.size() == tria.n_vertices());
00886 string strFileName = getDataFileName();
00887 m_List.push_back(strFileName);
00888 std::ofstream fout(strFileName.c_str());
00889 unsigned i;
00890 fout << "#Plotted For Vertices Values" << endl;
00891
00892 Triangulation<2>::active_cell_iterator cell = tria.begin_active(),
00893 endc = tria.end();
00894 for (; cell!=endc; ++cell)
00895 {
00896 for (i =0;i<GeometryInfo<2>::vertices_per_cell;i++)
00897 {
00898 for (unsigned j=0;j<2;j++)
00899 fout << cell->vertex(i)[j] << " ";
00900 fout << vSol(cell->vertex_index(i)) << "\n";
00901 }
00902 i=0;
00903 for (unsigned j=0;j<2;j++)
00904 fout << cell->vertex(i)[j] << " ";
00905 fout << vSol(cell->vertex_index(i)) << "\n";
00906
00907
00908 fout << "\n\n";
00909 }
00910 fout.close();
00911 plotPoints3d(strFileName,name,strComp);
00912 return;
00913
00914 }
00915
00916
00917 void GnuPlotDeal::plotDealSolutionAtVertices(Triangulation<2> &tria,const VecDouble &vSol,const VecDouble& vU0,const VecDouble &vU1,std::string name,std::string strComp,bool bScattered)
00918 {
00919 assert(vU0.size() == tria.n_vertices());
00920 assert(vU1.size() == tria.n_vertices());
00921 assert(vSol.size() == tria.n_vertices());
00922 unsigned i;
00923 string strFileName = getDataFileName();
00924 m_List.push_back(strFileName);
00925 std::ofstream fout(strFileName.c_str());
00926
00927 fout << "#Plotted For Vertices Values With Deformation" << endl;
00928
00929 Triangulation<2>::active_cell_iterator cell = tria.begin_active(),
00930 endc = tria.end();
00931 for (; cell!=endc; ++cell)
00932 {
00933 for (i =0;i<GeometryInfo<2>::vertices_per_cell;i++)
00934 {
00935 fout << cell->vertex(i)[0] + vU0(cell->vertex_index(i)) << " ";
00936 fout << cell->vertex(i)[1] + vU1(cell->vertex_index(i)) << " ";
00937 fout << vSol(cell->vertex_index(i)) << "\n";
00938 }
00939 if (!bScattered)
00940 {
00941 i=0;
00942 fout << cell->vertex(i)[0] + vU0(cell->vertex_index(i)) << " ";
00943 fout << cell->vertex(i)[1] + vU1(cell->vertex_index(i)) << " ";
00944 fout << vSol(cell->vertex_index(i)) << "\n\n\n";
00945 }
00946 }
00947 fout.close();
00948
00949 plotPoints3d(strFileName,name,strComp);
00950 return;
00951
00952 }
00953
00954
00955
00956
00957
00958
00959
00960
00961
00967 void GnuPlotDeal::plotDealCentralPointsOfFlaggedCells(const Triangulation<2> &tria,const VecBool &cVFlags,std::string name)
00968 {
00969 assert(cVFlags.size() == tria.n_active_cells());
00970 string strFileName = getDataFileName();
00971 m_List.push_back(strFileName);
00972 std::ofstream fout(strFileName.c_str());
00973
00974 fout << "#Plotted Central Points of Flagged Cells" << endl;
00975
00976 Triangulation<2>::active_cell_iterator cell = tria.begin_active(),
00977 endc = tria.end();
00978 bool bFileEmpty=true;
00979 for (; cell!=endc; ++cell)
00980 {
00981 if (cVFlags[cell->index()])
00982 {
00983 bFileEmpty=false;
00984 for (unsigned j=0;j<2;j++)
00985 fout << cell->barycenter()[j] << " ";
00986 fout << 0 << "\n";
00987 }
00988 }
00989 if (!bFileEmpty)
00990 {
00991 fout.close();
00992 if (2 == 2)
00993 plotPoints3d(strFileName,name,"with points");
00994 else
00995 plotPoints(strFileName,name,"with points");
00996 }
00997 return;
00998
00999
01000 }
01001
01002
01008 void GnuPlotDeal::plotDealFlaggedVertices(const Triangulation<2> &tria,const VecBool &vVFlags,std::string name,std::string strComp)
01009 {
01010 assert(vVFlags.size() == tria.n_vertices());
01011 string strFileName = getDataFileName();
01012 m_List.push_back(strFileName);
01013 std::ofstream fout(strFileName.c_str());
01014 bool empty=true;
01015
01016 fout << "#Plotted Vertice Points of Flagged Cells" << endl;
01017 Triangulation<2>::active_cell_iterator cell = tria.begin_active(),
01018 endc = tria.end();
01019 for (; cell!=endc; ++cell)
01020 {
01021 for (unsigned i =0;i<GeometryInfo<2>::vertices_per_cell;i++)
01022 {
01023 if (vVFlags[cell->vertex_index(i)])
01024 {
01025 for (unsigned j=0;j<2;j++)
01026 fout << cell->vertex(i)[j] << " ";
01027 fout << "0 \n";
01028 empty=false;
01029
01030 }
01031 }
01032 }
01033 fout << endl;
01034 fout.close();
01035 if (empty)
01036 return;
01037 if (2 == 2)
01038 plotPoints3d(strFileName,name,strComp);
01039 else
01040 plotPoints(strFileName,name,strComp);
01041 return;
01042
01043
01044
01045 }
01046
01047 void GnuPlotDeal::plotFunction(const Triangulation<2> &tria,Function2D &f,std::string name,std::string strComp)
01048 {
01049 string strFileName = getDataFileName();
01050 m_List.push_back(strFileName);
01051 std::ofstream fout(strFileName.c_str());
01052 Triangulation<2>::active_cell_iterator cell = tria.begin_active(),
01053 endc = tria.end();
01054 unsigned i;
01055 for (; cell!=endc; ++cell)
01056 {
01057 for (i =0;i<GeometryInfo<2>::vertices_per_cell;i++)
01058 {
01059 for (unsigned j=0;j<2;j++)
01060 fout << cell->vertex(i)[j] << " ";
01061 fout << f(cell->vertex(i)) << "\n";
01062 }
01063 i=0;
01064
01065 for (unsigned j=0;j<2;j++)
01066 fout << cell->vertex(i)[j] << " ";
01067 fout << f(cell->vertex(i)) << "\n";
01068
01069
01070 fout << "\n\n";
01071 }
01072 fout.close();
01073
01074
01075 if (2 == 1)
01076 {
01077 plotPoints(strFileName,name,strComp);
01078 }
01079 if (2 == 2)
01080 {
01081 plotPoints3d(strFileName,name,strComp);
01082 }
01083
01084 return;
01085
01086 }
01087
01088
01096 void GnuPlotDeal::plotDealSolution(DoFHandler<2> &dofs,const VecDouble &sol,const VecDouble vU0,const VecDouble vU1,std::string name,std::string strComp,unsigned deg)
01097 {
01098 assert(vU0.size() == dofs.get_tria().n_vertices());
01099 assert(vU1.size() == dofs.get_tria().n_vertices());
01100
01101 unsigned i=0;
01102 string strFileName = getDataFileName();
01103 m_List.push_back(strFileName);
01104 std::ofstream fout(strFileName.c_str());
01105 fout << "#Plotando x y f(x,y) no grau de liberdade " << deg << std::endl;
01106
01107 const unsigned int dofs_per_cell = dofs.get_fe().dofs_per_cell;
01108 std::vector<unsigned int> local_dof_indices (dofs_per_cell);
01109
01110 DoFHandler<2>::active_cell_iterator cell = dofs.begin_active(),
01111 endc = dofs.end();
01112
01113 for (; cell!=endc; ++cell)
01114 {
01115 cell->get_dof_indices (local_dof_indices);
01116 for (i =0;i<GeometryInfo<2>::vertices_per_cell;i++)
01117 {
01118 fout << cell->vertex(i)[0] + vU0(cell->vertex_index(i)) << " ";
01119 fout << cell->vertex(i)[1] + vU1(cell->vertex_index(i)) << " ";
01120 fout << sol(cell->vertex_dof_index(i,deg)) << "\n";
01121 }
01122 i=0;
01123 fout << cell->vertex(i)[0] + vU0(cell->vertex_index(i)) << " ";
01124 fout << cell->vertex(i)[1] + vU1(cell->vertex_index(i)) << " ";
01125 fout << sol(cell->vertex_dof_index(i,deg)) << "\n\n\n";
01126 }
01127 fout.close();
01128
01129 plotPoints3d(strFileName,name,strComp);
01130
01131 return;
01132
01133 }
01134
01135 void GnuPlotDeal::plotDealSolutionAtY(DoFHandler<2> &dofs,const VecDouble &sol,std::string name,std::string strComp,unsigned deg,double dAtY)
01136 {
01137
01138 string strFileName = getDataFileName();
01139 m_List.push_back(strFileName);
01140 std::ofstream fout(strFileName.c_str());
01141 fout << "#Plotando x y f(x,y) no grau de liberdade " << deg << std::endl;
01142
01143
01144 DoFHandler<2>::active_cell_iterator cell = dofs.begin_active(),
01145 endc = dofs.end();
01146
01147 double Tolerance = fabs(cell->vertex(BL_VERTEX)[1] - cell->vertex(UL_VERTEX)[1])/2.0;
01148 for (; cell!=endc; ++cell)
01149 {
01150 if ( fabs(cell->vertex(BL_VERTEX)[1]-dAtY) < Tolerance && fabs(cell->vertex(BR_VERTEX)[1]-dAtY) < Tolerance)
01151 {
01152 fout << cell->vertex(BL_VERTEX)[0] << " " << sol(cell->vertex_dof_index(BL_VERTEX,deg)) << endl;
01153 fout << cell->vertex(BR_VERTEX)[0] << " " << sol(cell->vertex_dof_index(BR_VERTEX,deg)) << endl << endl << endl;
01154 }
01155 else if ( fabs(cell->vertex(UL_VERTEX)[1]-dAtY) < Tolerance && fabs(cell->vertex(UR_VERTEX)[1]-dAtY) < Tolerance)
01156 {
01157 fout << cell->vertex(UL_VERTEX)[0] << " " << sol(cell->vertex_dof_index(UL_VERTEX,deg)) << endl;
01158 fout << cell->vertex(UR_VERTEX)[0] << " " << sol(cell->vertex_dof_index(UR_VERTEX,deg)) << endl << endl << endl;
01159 }
01160
01161 }
01162 fout.close();
01163
01164 plotPoints(strFileName,name,strComp);
01165
01166 return;
01167
01168 }
01169
01170 void GnuPlotDeal::plotFunction1D(Function1D &f,double dStart,double dEnd,int nPoints,std::string name,std::string strComp,unsigned deg)
01171 {
01172 double dP;
01173 double dx = (dEnd - dStart)/nPoints;
01174
01175
01176 string strFileName = getDataFileName();
01177 m_List.push_back(strFileName);
01178 std::ofstream fout(strFileName.c_str());
01179 fout << "#Plotting Function1D" << std::endl;
01180
01181 for (dP = dStart;dP<=dEnd;dP+=dx)
01182 {
01183 Point1D P(dP);
01184 fout << dP << " " << f(P,deg) << endl;
01185 }
01186 fout.close();
01187 plotPoints(strFileName,name,strComp);
01188 return;
01189
01190 }
01191
01192 void GnuPlotDeal::plotDealSolutionAtVerticesAtY(Triangulation<2> &tria,const VecDouble &vSol,double dAtY,std::string name,std::string strComp)
01193 {
01194 string strFileName = getDataFileName();
01195 m_List.push_back(strFileName);
01196 std::ofstream fout(strFileName.c_str());
01197 fout << "#Plotando x y f(x,y) nos vertices com Y == " << dAtY << std::endl;
01198
01199
01200 Triangulation<2>::active_cell_iterator cell = tria.begin_active(),
01201 endc = tria.end();
01202
01203 double Tolerance = fabs(cell->vertex(BL_VERTEX)[1] - cell->vertex(UL_VERTEX)[1])/2.0;
01204 for (; cell!=endc; ++cell)
01205 {
01206 if ( fabs(cell->vertex(BL_VERTEX)[1]-dAtY) < Tolerance && fabs(cell->vertex(BR_VERTEX)[1]-dAtY) < Tolerance)
01207 {
01208 fout << cell->vertex(BL_VERTEX)[0] << " " << vSol(cell->vertex_index(BL_VERTEX)) << endl;
01209 fout << cell->vertex(BR_VERTEX)[0] << " " << vSol(cell->vertex_index(BR_VERTEX)) << endl << endl << endl;
01210 }
01211 else if ( fabs(cell->vertex(UL_VERTEX)[1]-dAtY) < Tolerance && fabs(cell->vertex(UR_VERTEX)[1]-dAtY) < Tolerance)
01212 {
01213 fout << cell->vertex(UL_VERTEX)[0] << " " << vSol(cell->vertex_index(UL_VERTEX)) << endl;
01214 fout << cell->vertex(UR_VERTEX)[0] << " " << vSol(cell->vertex_index(UR_VERTEX)) << endl << endl << endl;
01215 }
01216
01217 }
01218 fout.close();
01219
01220 plotPoints(strFileName,name,strComp);
01221
01222 return;
01223
01224 }
01225
01226
01233 void GnuPlotDeal::plotDealSolutionAtX(DoFHandler<2> &dofs,const VecDouble &sol,std::string name,std::string strComp,unsigned deg,double dAtX)
01234 {
01235 string strFileName = getDataFileName();
01236 m_List.push_back(strFileName);
01237 std::ofstream fout(strFileName.c_str());
01238 fout << "#Plotando x y f(x,y) no grau de liberdade " << deg << std::endl;
01239
01240
01241 DoFHandler<2>::active_cell_iterator cell = dofs.begin_active(),
01242 endc = dofs.end();
01243
01244 double Tolerance = fabs(cell->vertex(BL_VERTEX)[1] - cell->vertex(UL_VERTEX)[1])/2.0;
01245 for (; cell!=endc; ++cell)
01246 {
01247 if ( fabs(cell->vertex(BL_VERTEX)[0]-dAtX) < Tolerance && fabs(cell->vertex(UL_VERTEX)[0]-dAtX) < Tolerance)
01248 {
01249 fout << cell->vertex(BL_VERTEX)[1] << " " << sol(cell->vertex_dof_index(BL_VERTEX,deg)) << endl;
01250 fout << cell->vertex(UL_VERTEX)[1] << " " << sol(cell->vertex_dof_index(UL_VERTEX,deg)) << endl << endl << endl;
01251 }
01252 else if ( fabs(cell->vertex(BR_VERTEX)[0]-dAtX) < Tolerance && fabs(cell->vertex(UR_VERTEX)[0]-dAtX) < Tolerance)
01253 {
01254 fout << cell->vertex(BR_VERTEX)[1] << " " << sol(cell->vertex_dof_index(BR_VERTEX,deg)) << endl;
01255 fout << cell->vertex(UR_VERTEX)[1] << " " << sol(cell->vertex_dof_index(UR_VERTEX,deg)) << endl << endl << endl;
01256 }
01257
01258 }
01259 fout.close();
01260
01261 plotPoints(strFileName,name,strComp);
01262
01263 return;
01264
01265 }
01266
01267
01268
01269
01270
01271
01272
01273
01274 void GnuPlotDeal::plotShapeFunctionsAtCell(const FiniteElement<2> &fe,DoFHandler<2>::cell_iterator &cell,unsigned nPoints)
01275 {
01276 double dx = 1.0/nPoints;
01277 char str[500];
01278 Point<2> p(0,0);
01279
01280 int n_shape_functions = fe.n_dofs_per_cell();
01281 std::string strFileName;
01282 std::ofstream fout;
01283
01284
01285 p[0]=p[1]=0;
01286 std::vector<Point<2> > lstPts;
01287 unsigned isoLinesCount;
01288 for (unsigned i=0;i<=nPoints;i++,p[0]+=dx)
01289 {
01290 for (unsigned j=0;j<=nPoints;j++,p[1]+=dx)
01291 lstPts.push_back(p);
01292 p[1]=0;
01293 }
01294 Quadrature<2> quadrature(lstPts);
01295 FEValues<2> fe_values (fe, quadrature,
01296 UpdateFlags(update_values|update_q_points));
01297
01298 fe_values.reinit(cell);
01299 for (int iShape=0;iShape<n_shape_functions;iShape++)
01300 {
01301 for (unsigned component=0;component<fe.get_nonzero_components(iShape).size();component++)
01302 {
01303
01304 if (fe.get_nonzero_components(iShape)[component])
01305 {
01306 sprintf(str,"Shape: %d Comp: %d",iShape,component);
01307 strFileName=newDataFile(fout);
01308 isoLinesCount = 1;
01309 for (unsigned qPoint=0;qPoint < fe_values.n_quadrature_points;qPoint++,isoLinesCount++)
01310 {
01311 fout << fe_values.quadrature_point(qPoint)[0] << " " << fe_values.quadrature_point(qPoint)[1] << " " <<fe_values.shape_value_component(iShape,qPoint,component) << "\n";
01312 if (isoLinesCount > nPoints)
01313 {
01314 isoLinesCount = 0;
01315 fout << "\n";
01316 }
01317 }
01318
01319 fout.close();
01320 plotPoints3d(strFileName,str,"with lines");
01321 }
01322 }
01323 }
01324 }
01325
01326
01334 void GnuPlotDeal::plotDealGradComponent(DoFHandler<2> &dofs,const VecDouble &sol,unsigned gradComp,std::string name,std::string strComp,VecDouble cVWeights)
01335 {
01336 assert(cVWeights.size() == 0 || cVWeights.size() == dofs.get_tria().n_active_cells());
01337
01338 string strFileName = getDataFileName();
01339 m_List.push_back(strFileName);
01340 std::ofstream fout(strFileName.c_str());
01341 double factor;
01342 const unsigned int dofs_per_cell = dofs.get_fe().dofs_per_cell;
01343 std::vector<unsigned int> local_dof_indices (dofs_per_cell);
01344
01345 DoFHandler<2>::active_cell_iterator cell = dofs.begin_active(),
01346 endc = dofs.end();
01347
01348 std::vector< Point<2> > points;
01349 points.push_back(Point<2>(0,0));
01350 points.push_back(Point<2>(1,0));
01351 points.push_back(Point<2>(1,1));
01352 points.push_back(Point<2>(0,1));
01353 points.push_back(Point<2>(0,0));
01354
01355
01356 Quadrature<2> quad(points);
01357
01358 FEValues<2> fe_values (dofs.get_fe(), quad,
01359 UpdateFlags(update_gradients|update_q_points));
01360 for (; cell!=endc; ++cell)
01361 {
01362 cell->get_dof_indices (local_dof_indices);
01363 fe_values.reinit(cell);
01364 if (cVWeights.size() != 0)
01365 factor = cVWeights(cell->index());
01366 else
01367 factor = 1;
01368 for (unsigned qPoint=0;qPoint<quad.n_quadrature_points;qPoint++)
01369 {
01370 double result=0;
01371 for (unsigned i=0;i<dofs_per_cell;i++)
01372 result += factor*sol(local_dof_indices[i])*fe_values.shape_grad(i,qPoint)[gradComp];
01373 fout << fe_values.quadrature_point(qPoint)[0] << " " << fe_values.quadrature_point(qPoint)[1] << " " << result << endl;
01374 }
01375 fout << endl << endl;
01376 }
01377 fout.close();
01378 plotPoints3d(strFileName,name,strComp);
01379
01380 }
01381
01382
01387 void GnuPlotDeal::plotDealSolutionComponent(DoFHandler<2> &dofs,const VecDouble &sol,unsigned comp,std::string name,std::string strComp)
01388 {
01389 string strFileName = getDataFileName();
01390 m_List.push_back(strFileName);
01391 std::ofstream fout(strFileName.c_str());
01392
01393 const unsigned int dofs_per_cell = dofs.get_fe().dofs_per_cell;
01394 std::vector<unsigned int> local_dof_indices (dofs_per_cell);
01395
01396 DoFHandler<2>::active_cell_iterator cell = dofs.begin_active(),
01397 endc = dofs.end();
01398
01399 std::vector< Point<2> > points;
01400 points.push_back(Point<2>(0,0));
01401 points.push_back(Point<2>(1,0));
01402 points.push_back(Point<2>(1,1));
01403 points.push_back(Point<2>(0,1));
01404 points.push_back(Point<2>(0,0));
01405
01406 Quadrature<2> quad(points);
01407
01408 FEValues<2> fe_values (dofs.get_fe(), quad,
01409 UpdateFlags(update_values|update_q_points));
01410 for (; cell!=endc; ++cell)
01411 {
01412 cell->get_dof_indices (local_dof_indices);
01413 fe_values.reinit(cell);
01414 for (unsigned qPoint=0;qPoint<quad.n_quadrature_points;qPoint++)
01415 {
01416 double result=0;
01417 for (unsigned i=0;i<dofs_per_cell;i++)
01418 result += sol(local_dof_indices[i])*fe_values.shape_value_component(i,qPoint,comp);
01419 fout << fe_values.quadrature_point(qPoint)[0] << " " << fe_values.quadrature_point(qPoint)[1] << " " << result << endl;
01420 }
01421 fout << endl << endl;
01422 }
01423 fout.close();
01424 plotPoints3d(strFileName,name,strComp);
01425
01426 }
01427
01428
01432 void GnuPlotDeal::plotMapInVectors(const VecDouble &v1,const VecDouble &v2,std::string name,std::string strComp)
01433 {
01434 if (v1.size()!=v2.size()) abort();
01435
01436 double dd[v1.size()*2];
01437 unsigned j=0;
01438 for (unsigned i=0;i<v1.size();i++)
01439 {
01440 assert(j<v1.size()*2);
01441 dd[j++]=v1(i);
01442 dd[j++]=v2(i);
01443 }
01444 this->plotPoints(dd,v1.size(),2,name,"",strComp,true,false);
01445 }
01446
01447
01448
01456 void GnuPlotDeal::plotDealSolutionAtVerticesAtX(Triangulation<2> &tria,const VecDouble &vSol,double dAtX,std::string name,std::string strComp)
01457 {
01458 string strFileName = getDataFileName();
01459 m_List.push_back(strFileName);
01460 std::ofstream fout(strFileName.c_str());
01461 fout << "#Plotando x y f(x,y) em " << dAtX << std::endl;
01462
01463 Triangulation<2>::active_cell_iterator cell = tria.begin_active();
01464 Triangulation<2>::active_cell_iterator endc = tria.end();
01465
01466 double Tolerance = fabs(cell->vertex(BL_VERTEX)[1] - cell->vertex(UL_VERTEX)[1])/2.0;
01467 for (; cell!=endc; ++cell)
01468 {
01469 if ( fabs(cell->vertex(BL_VERTEX)[0]-dAtX) < Tolerance && fabs(cell->vertex(UL_VERTEX)[0]-dAtX) < Tolerance)
01470 {
01471 fout << cell->vertex(BL_VERTEX)[1] << " " << vSol(cell->vertex_index(BL_VERTEX)) << endl;
01472 fout << cell->vertex(UL_VERTEX)[1] << " " << vSol(cell->vertex_index(UL_VERTEX)) << endl << endl << endl;
01473 }
01474 else if ( fabs(cell->vertex(BR_VERTEX)[0]-dAtX) < Tolerance && fabs(cell->vertex(UR_VERTEX)[0]-dAtX) < Tolerance)
01475 {
01476 fout << cell->vertex(BR_VERTEX)[1] << " " << vSol(cell->vertex_index(BR_VERTEX)) << endl;
01477 fout << cell->vertex(UR_VERTEX)[1] << " " << vSol(cell->vertex_index(UR_VERTEX)) << endl << endl << endl;
01478 }
01479
01480 }
01481 fout.close();
01482
01483 plotPoints(strFileName,name,strComp);
01484
01485 return;
01486
01487 }
01488
01489
01490
01491
01499 void GnuPlotDeal::plotDealSolutionCentralPointsAtY(Triangulation<2> &tria,const VecDouble &vSol,double dAtY,std::string name,std::string strComp)
01500 {
01501 assert(vSol.size() == tria.n_active_cells());
01502 string strFileName = getDataFileName();
01503 m_List.push_back(strFileName);
01504 std::ofstream fout(strFileName.c_str());
01505 fout << "#Plotando x y f(x,y) nos vertices com Y == " << dAtY << std::endl;
01506
01507
01508 Triangulation<2>::active_cell_iterator cell = tria.begin_active(),
01509 endc = tria.end(),cellFront;
01510
01511 double Tolerance = fabs(cell->vertex(BL_VERTEX)[1] - cell->vertex(UL_VERTEX)[1])*0.7;
01512 for (; cell!=endc; ++cell)
01513 {
01514 if ( fabs(cell->vertex(BL_VERTEX)[1]-dAtY) < Tolerance && fabs(cell->vertex(BR_VERTEX)[1]-dAtY) < Tolerance)
01515 {
01516 if (cell->face(LEFT_CELL)->at_boundary())
01517 break;
01518 }
01519 }
01520 assert(cell->face(LEFT_CELL)->at_boundary());
01521
01522 while(! cell->face(RIGHT_CELL)->at_boundary())
01523 {
01524 cellFront = cell->neighbor(RIGHT_CELL);
01525 fout << cell->vertex(BR_VERTEX)[0] << " " << ( vSol(cell->index()) + vSol(cellFront->index()) )/2.0 << "\n";
01526 cell=cellFront;
01527 }
01528 fout.close();
01529
01530 plotPoints(strFileName,name,strComp);
01531
01532 return;
01533
01534
01535 }
01536
01544 void GnuPlotDeal::plotCentralValuesAtFixedZ(Triangulation<3> &tria, const VecDouble &cSol,string name, string strComp,double atZ)
01545 {
01546 assert(cSol.size() == tria.n_active_cells());
01547 string strFileName = getDataFileName();
01548 m_List.push_back(strFileName);
01549 std::ofstream fout(strFileName.c_str());
01550
01551 fout << "#Central Values at Z = " << atZ << endl;
01552
01553 Triangulation<3>::active_cell_iterator cell = tria.begin_active(),
01554 endc = tria.end();
01555 for (; cell!=endc; ++cell)
01556 {
01557 if ( atZ > cell->vertex(VERTEX_000)[2] && atZ < cell->vertex(VERTEX_001)[2])
01558 {
01559 fout << cell->vertex(VERTEX_000)[0] << " " << cell->vertex(VERTEX_000)[1] << " " << cSol(cell->index()) << "\n";
01560 fout << cell->vertex(VERTEX_100)[0] << " " << cell->vertex(VERTEX_100)[1] << " " << cSol(cell->index()) << "\n";
01561 fout << cell->vertex(VERTEX_110)[0] << " " << cell->vertex(VERTEX_110)[1] << " " << cSol(cell->index()) << "\n";
01562 fout << cell->vertex(VERTEX_010)[0] << " " << cell->vertex(VERTEX_010)[1] << " " << cSol(cell->index()) << "\n";
01563 fout << cell->vertex(VERTEX_000)[0] << " " << cell->vertex(VERTEX_000)[1] << " " << cSol(cell->index()) << "\n";
01564 fout << endl << endl;
01565 }
01566 }
01567 fout.close();
01568
01569 plotPoints3d(strFileName,name,strComp);
01570
01571 return;
01572
01573
01574
01575 }
01576
01577
01585 void GnuPlotDeal::plotVerticeValuesAtFixedZ(Triangulation<3> &tria, const VecDouble &vSol,string name, string strComp,double atZ)
01586 {
01587 assert(vSol.size() == tria.n_vertices());
01588 string strFileName = getDataFileName();
01589 m_List.push_back(strFileName);
01590 std::ofstream fout(strFileName.c_str());
01591
01592 fout << "#Central Values at Z = " << atZ << endl;
01593
01594 Triangulation<3>::active_cell_iterator cell = tria.begin_active(),
01595 endc = tria.end();
01596 for (; cell!=endc; ++cell)
01597 {
01598 if ( atZ > cell->vertex(VERTEX_000)[2] && atZ < cell->vertex(VERTEX_001)[2])
01599 {
01600 fout << cell->vertex(VERTEX_000)[0] << " " << cell->vertex(VERTEX_000)[1] << " " << vSol(cell->vertex_index(VERTEX_000)) << "\n";
01601 fout << cell->vertex(VERTEX_100)[0] << " " << cell->vertex(VERTEX_100)[1] << " " << vSol(cell->vertex_index(VERTEX_100)) << "\n";
01602 fout << cell->vertex(VERTEX_110)[0] << " " << cell->vertex(VERTEX_110)[1] << " " << vSol(cell->vertex_index(VERTEX_110)) << "\n";
01603 fout << cell->vertex(VERTEX_010)[0] << " " << cell->vertex(VERTEX_010)[1] << " " << vSol(cell->vertex_index(VERTEX_010)) << "\n";
01604 fout << cell->vertex(VERTEX_000)[0] << " " << cell->vertex(VERTEX_000)[1] << " " << vSol(cell->vertex_index(VERTEX_000)) << "\n";
01605 fout << endl << endl;
01606 }
01607 }
01608 fout.close();
01609
01610 plotPoints3d(strFileName,name,strComp);
01611
01612 return;
01613
01614 }
01615
01616
01617
01618 void GnuPlotDeal::plotVerticeValuesAtFixedZ(OrthoMesh &mesh, const VecDouble &vSol,string name, string strComp,double atZ)
01619 {
01620 assert(vSol.size() == mesh.numVertices());
01621 string strFileName = getDataFileName();
01622 m_List.push_back(strFileName);
01623 std::ofstream fout(strFileName.c_str());
01624
01625 fout << "#Central Values at Z = " << atZ << endl;
01626
01627 OrthoMesh::Cell_It cell = mesh.begin_cell(),
01628 endc = mesh.end_cell();
01629 for (; cell!=endc; cell++)
01630 {
01631 if ( atZ >= cell->vertex(VERTEX_000)[2] && atZ < cell->vertex(VERTEX_001)[2])
01632 {
01633 fout << cell->vertex(VERTEX_000)[0] << " " << cell->vertex(VERTEX_000)[1] << " " << vSol(cell->vertex_index(VERTEX_000)) << "\n";
01634 fout << cell->vertex(VERTEX_100)[0] << " " << cell->vertex(VERTEX_100)[1] << " " << vSol(cell->vertex_index(VERTEX_100)) << "\n";
01635 fout << cell->vertex(VERTEX_110)[0] << " " << cell->vertex(VERTEX_110)[1] << " " << vSol(cell->vertex_index(VERTEX_110)) << "\n";
01636 fout << cell->vertex(VERTEX_010)[0] << " " << cell->vertex(VERTEX_010)[1] << " " << vSol(cell->vertex_index(VERTEX_010)) << "\n";
01637 fout << cell->vertex(VERTEX_000)[0] << " " << cell->vertex(VERTEX_000)[1] << " " << vSol(cell->vertex_index(VERTEX_000)) << "\n";
01638 fout << endl << endl;
01639 }
01640 }
01641 fout.close();
01642
01643 plotPoints3d(strFileName,name,strComp);
01644
01645 return;
01646
01647 }
01648
01649
01650
01651 void GnuPlotDeal::plotLineAlongZ(OrthoMesh &mesh,const VecDouble &vSol,string name, string strComp,double X,double Y)
01652 {
01653 assert(vSol.size() == mesh.numVertices());
01654 string strFileName = getDataFileName();
01655 m_List.push_back(strFileName);
01656 std::ofstream fout(strFileName.c_str());
01657
01658 fout << "#Values along Z (" << X << ", " << Y << ", Z)\n" ;
01659
01660 OrthoMesh::Cell_It cell = mesh.begin_cell(),
01661 endc = mesh.end_cell();
01662 bool found=false;
01663 for (; cell!=endc; cell++)
01664 {
01665 if ( X >= cell->vertex(VERTEX_000)[0] && X <= cell->vertex(VERTEX_110)[0] &&
01666 Y >= cell->vertex(VERTEX_000)[1] && Y <= cell->vertex(VERTEX_110)[1])
01667 {
01668 found=true;
01669 break;
01670 }
01671 }
01672 if (!found)
01673 {
01674 printf("GnuPlot Warning: Point %g %g does not belong to domain\n",X,Y);
01675 fout.close();
01676 return;
01677 }
01678
01679 fout << cell->vertex(VERTEX_000)[2] << " " << vSol(cell->vertex_index(VERTEX_000)) << std::endl;
01680 fout << cell->vertex(VERTEX_001)[2] << " " << vSol(cell->vertex_index(VERTEX_001)) << std::endl;
01681 while (cell->advance(UP_CELL))
01682 {
01683 fout << cell->vertex(VERTEX_001)[2] << " " << vSol(cell->vertex_index(VERTEX_001)) << std::endl;
01684 }
01685
01686 plotPoints(strFileName,name,strComp);
01687
01688 return;
01689
01690 }
01691
01692
01693 void GnuPlotDeal::plotLineAlongX(OrthoMesh &mesh,const VecDouble &vSol,string name, string strComp,double Y,double Z)
01694 {
01695 assert(vSol.size() == mesh.numVertices());
01696 string strFileName = getDataFileName();
01697 m_List.push_back(strFileName);
01698 std::ofstream fout(strFileName.c_str());
01699
01700
01701 writeDataAlongX(mesh,vSol,Y,Z,fout);
01702 plotPoints(strFileName,name,strComp);
01703
01704 return;
01705
01706 }
01707
01708
01709 void GnuPlotDeal::writeDataAlongX(OrthoMesh &mesh,const VecDouble &vSol,double Y,double Z,std::ostream &fout)
01710 {
01711 assert(vSol.size() == mesh.numVertices());
01712 fout << "#Values along X (" << "X" << ", " << Y << ", " << Z << " )\n" ;
01713
01714 OrthoMesh::Cell_It cell = mesh.begin_cell(),
01715 endc = mesh.end_cell();
01716 bool found=false;
01717 for (; cell!=endc; cell++)
01718 {
01719 if ( Y >= cell->vertex(VERTEX_000)[1] && Y <= cell->vertex(VERTEX_111)[1] &&
01720 Z >= cell->vertex(VERTEX_000)[2] && Y <= cell->vertex(VERTEX_111)[2])
01721 {
01722 found=true;
01723 break;
01724 }
01725 }
01726 if (!found)
01727 {
01728 printf("GnuPlot Warning: Point %g %g does not belong to domain\n",Y,Z);
01729 fout.flush();
01730 return;
01731 }
01732
01733 fout << cell->vertex(VERTEX_000)[0] << " " << vSol(cell->vertex_index(VERTEX_000)) << std::endl;
01734 fout << cell->vertex(VERTEX_100)[0] << " " << vSol(cell->vertex_index(VERTEX_100)) << std::endl;
01735 while (cell->advance(RIGHT_CELL))
01736 {
01737 fout << cell->vertex(VERTEX_100)[0] << " " << vSol(cell->vertex_index(VERTEX_100)) << std::endl;
01738 }
01739
01740 return;
01741
01742
01743 }
01744
01745
01746 void GnuPlotDeal::writeDataAlongZ(OrthoMesh &mesh,const VecDouble &vSol,double X,double Y,std::ostream &fout)
01747 {
01748 assert(vSol.size() == mesh.numVertices());
01749
01750 fout << "#Values along Z (" << X << ", " << Y << ", " << "Z" << " )\n" ;
01751
01752 OrthoMesh::Cell_It cell = mesh.begin_cell(),
01753 endc = mesh.end_cell();
01754 bool found=false;
01755 for (; cell!=endc; cell++)
01756 {
01757 if ( X >= cell->vertex(VERTEX_000)[0] && X <= cell->vertex(VERTEX_111)[0] &&
01758 Y >= cell->vertex(VERTEX_000)[1] && Y <= cell->vertex(VERTEX_111)[1])
01759 {
01760 found=true;
01761 break;
01762 }
01763 }
01764 if (!found)
01765 {
01766 printf("GnuPlot Warning: Point %g %g does not belong to domain\n",X,Y);
01767 fout.flush();
01768 return;
01769 }
01770
01771 fout << cell->vertex(VERTEX_000)[2] << " " << vSol(cell->vertex_index(VERTEX_000)) << std::endl;
01772 fout << cell->vertex(VERTEX_001)[2] << " " << vSol(cell->vertex_index(VERTEX_001)) << std::endl;
01773 while (cell->advance(UP_CELL))
01774 {
01775 fout << cell->vertex(VERTEX_001)[2] << " " << vSol(cell->vertex_index(VERTEX_001)) << std::endl;
01776 }
01777
01778 return;
01779
01780
01781
01782 }
