Mercurial > hg > kwantix

changeset 24:ad9e3d28ce9b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Implement op<< for bvp::interpolator
author Jordi Gutiérrez Hermoso <jordigh@gmail.com>
date Fri, 21 Aug 2009 16:19:58 -0500
parents acefa48d4b4d
children 6c6003bcad16
files Makefile data/linear_wave.m include/interpolator.hpp interpolator.cpp linalg.cpp main.cpp results/plotresults_vtk.m vtk/terrain_test.cpp vtk/terrain_test.tcl
diffstat 9 files changed, 298 insertions(+), 25 deletions(-) [+]
line wrap: on
line diff
--- a/Makefile
+++ b/Makefile
@@ -1,6 +1,6 @@
 CPP = g++
-LINKING = -lgsl -lgslcblas 
-CFLAGS = -O3
+LINKING = -lgsl -lgslcblas
+CFLAGS = -pg -O2 -g
 OPTIONS = -Wall -pedantic  -W -Werror -Wconversion -Wshadow \
 	  -Wpointer-arith -Wcast-qual -Wcast-align -Wwrite-strings  \
 	  -fshort-enums -fno-common -Wfatal-errors 
--- a/data/linear_wave.m
+++ b/data/linear_wave.m
@@ -1,5 +1,5 @@
-n = 20
-m = 6;
+n = 19
+m = 8;
 r = linspace(0,1,n);
 
 rrintr = 0;
@@ -46,9 +46,9 @@
 ##nonradial
 h0 = [];
  for i = 1:columns(xxintr)
-   h0(i) = fr( (xxintr(i)-0.25)^2 + yyintr(i)^2, 0.025);
+   h0(i) = fr( (xxintr(i)-0.25)^2 + yyintr(i)^2, 0.040);
  endfor
- h0 = 4*h0(:)+0.1;
+ h0 = 4*h0(:);
 
 h0 = [h0];
 
--- a/include/interpolator.hpp
+++ b/include/interpolator.hpp
@@ -14,6 +14,7 @@
 #include <vector>
 #include <map>
 #include <boost/shared_ptr.hpp>
+#include <iostream>
 #include "bvp.hpp"
 #include "linalg.hpp"
 #include "func.hpp"
@@ -188,9 +189,7 @@
      *   interpolator is defined and the \f$n^{th}\f$ column is the
      *   values that the interpolator takes on that point.
      */
-    friend
-    template<typename RBF>
-    linalg::matrix& operator<< (ostream& os, 
+    void into_os(std::ostream& os) const;
     
   private:
     ///Once the matrix is defined, this function inverts it.
@@ -278,7 +277,15 @@
   {
     return u*a;
   }
-  
+
+  /// Stream insertion
+  template <typename RBF>
+  std::ostream& operator<< (std::ostream& os, const interpolator<RBF>& u)
+  {
+    u.into_os(os);
+    return os;
+  }
+
   
 }
 
--- a/interpolator.cpp
+++ b/interpolator.cpp
@@ -68,7 +68,8 @@
       else
       {
         badArgument exc;
-        exc.reason = "The interpolation data contains points not in the given domain.";
+        exc.reason = "The interpolation data contains points not in "
+                     "the given domain.";
         exc.line = __LINE__;
         exc.file = __FILE__;
         throw exc;
@@ -697,7 +698,50 @@
     return out;
   }
   
+  //****************** I/O *************************************
+  
+  template<typename RBF>
+  void interpolator<RBF>::into_os(std::ostream& os) const
+  {
+    linalg::vector values;
+    size_t i;
+    set<point>::const_iterator I;
+    
+    //If we already have precomputed the values...
+    std::vector<size_t> alpha;
+    if(precomp_values_vec.count(alpha)){
+      values = precomp_values_vec[alpha];
+    }
+    else{
+      linalg::vector tmp(n+m);
+      for(i = 1, I=thebvp->get_domain()->get_interior().begin() ; 
+          i<=n; i++, I++)
+      {
+        tmp(i) = at(*I);
+      }
 
+      for(i = n+1, I=thebvp->get_domain()->get_boundary().begin() ; 
+          i<=n+m; i++, I++)
+      {
+        tmp(i) = at(*I);
+      }
+
+      values = tmp;
+    }
+    
+    for(i = 1, I=thebvp->get_domain()->get_interior().begin() ; 
+        i<=n; i++, I++)
+    {
+      os << *I << " " << values(i) << std::endl;
+    }
+
+    for(i = n+1, I=thebvp->get_domain()->get_boundary().begin() ; 
+        i<=n+m; i++, I++)
+    {
+      os << *I << " " << values(i) << std::endl;
+    }
+ 
+  }
 
   //Instantiations
   using namespace rbf;
--- a/linalg.cpp
+++ b/linalg.cpp
@@ -770,9 +770,8 @@
     os.setf(std::ios::scientific);
     os << std::setprecision(int(v.precision()) );
     for(size_t i = 1; i <= v.size(); i++){
-      os << " " << std::setw(int(v.precision()+6) ) << v(i) << std::endl;
+      os << " " << std::setw(int(v.precision()+6) ) << v(i) ;
     }
-    os << std::endl;
      return os;
    }
 
--- a/main.cpp
+++ b/main.cpp
@@ -28,7 +28,7 @@
   wave_op(double dt_in): dt2(dt_in*dt_in){};
   double at(const realfunc& f, const point& p) const
   {
-    return dt2*L(f,p) + f(p);
+    return dt2*L(f,p) - f(p);
   }
 private:
   Laplacian L;
@@ -40,6 +40,9 @@
   double at(const point&) const{return 0;};
 };
 
+void save_iteration(const interpolator<RBF_TYPE>& u,
+                    size_t n);
+
 int main()
 {  
   gsl_set_error_handler(&error_handling::errorHandler);
@@ -52,8 +55,7 @@
       
     shared_ptr<domain> Omega(new domain("data/circ_intr.matrix",
                                         "data/circ_bdry.matrix",
-                                        "data/circ_nrml.matrix"));
-    
+                                        "data/circ_nrml.matrix"));    
     zero_func g;
 
     //timestep
@@ -70,9 +72,10 @@
     interpolator<RBF_TYPE> u0 = u, u1 = u;
 
     //Main loop
-    size_t maxiter = 1000;
+    size_t maxiter = 500;
     for(size_t n = 1; n <= maxiter; n++)
     {
+      save_iteration(u,n);
       u0 = u1;
       u1 = u;
       cout << "Iteration: " << n << endl;
@@ -87,4 +90,16 @@
 
 }
 
-
+void 
+save_iteration(const interpolator<RBF_TYPE>& u,
+               size_t n)
+{
+  using namespace std;
+  stringstream ss;
+  string filename = "results/u";
+  ss << filename << setw(5) << setfill('0') << n << ".map";
+  ss >> filename;
+  ofstream ofs(filename.c_str());
+  ofs << u;
+  ofs.close();
+}
--- a/results/plotresults_vtk.m
+++ b/results/plotresults_vtk.m
@@ -1,6 +1,6 @@
-skip = 5;
-for i = skip:skip:100000
-  name = "h";
+skip = 1;
+for i = skip:skip:500
+  name = "u";
   if(i < 10000)
     name = strcat(name, "0");
   endif
@@ -14,15 +14,16 @@
     name = strcat(name, "0");
   endif
   name = strcat(name,num2str(i));
-  vtk_xlim([-1,1]);
-  vtk_ylim([-1,1]);
-  vtk_zlim([0,0.2]);
+  #vtk_xlim([-1,1]);
+  #vtk_ylim([-1,1]);
+  #vtk_zlim([0,0.2]);
   filename = strcat(name,".map");
   load(filename);
   h = eval(name);
   T = delaunay(h(:,1),h(:,2));
   vtk_clear; vtk_trisurf(T,h(:,1),h(:,2),h(:,3));
-  ##vtk_trimesh(T,h(:,1),h(:,2),h(:,3));
+  ##foo = vtk_trimesh(T,h(:,1),h(:,2),h(:,3));
+  vtk_axis('off');
   vtk_print(strcat(name,".png"),"-dpng");
   clear eval(name);
 endfor
new file mode 100644
--- /dev/null
+++ b/vtk/terrain_test.cpp
@@ -0,0 +1,112 @@
+#include <vtkPoints.h>
+#include <vtkPointData.h>
+#include <vtkDelaunay2D.h>
+#include <vtkLookupTable.h>
+
+#include <vtkPolyData.h>
+#include <vtkPolyDataNormals.h>
+#include <vtkPolyDataMapper.h>
+
+#include <vtkWarpScalar.h>
+
+#include <vtkCamera.h>
+#include <vtkActor.h>
+#include <vtkRenderWindow.h>
+#include <vtkRenderWindowInteractor.h>
+#include <vtkRenderer.h>
+
+#include <vtkInteractorStyleTerrain.h>
+
+#include <vtkSmartPointer.h>
+#include <vtkDoubleArray.h>
+ 
+ 
+void TriangulateTerrain();
+
+int main (int argc, char ** argv)
+{
+	TriangulateTerrain();	
+	return 0;
+}
+  
+void TriangulateTerrain()
+{
+ 
+	vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
+  vtkSmartPointer<vtkDoubleArray> scalars = vtkSmartPointer<vtkDoubleArray>::New();
+ 
+	unsigned int GridSize = 50;
+	for(double x = -3; x < 3; x += 6.0/GridSize)
+	{
+		for(double y = -3; y < 3; y += 6.0/GridSize)
+		{
+      double z = sin(x*y);
+			vtkIdType idx = points->InsertNextPoint(x, y, z);
+      scalars -> InsertTuple(idx, &z);
+		}
+ 
+	}
+ 
+	//add the grid points to a polydata object
+	vtkPolyData* polydata = vtkPolyData::New();
+	polydata->SetPoints(points);
+  polydata->GetPointData() -> SetScalars(scalars);
+ 
+	//triangulate the grid points
+	vtkDelaunay2D* delaunay = vtkDelaunay2D::New();
+	delaunay->SetInput(polydata);
+	delaunay->Update();
+
+  //Normals for Gourad shading
+  vtkPolyDataNormals* normals = vtkPolyDataNormals::New();
+  normals -> SetInputConnection(delaunay ->GetOutputPort() );
+  normals -> SetFeatureAngle(60);
+
+  //Set the colours for the rendering
+  vtkLookupTable* lut = vtkLookupTable::New();
+  lut -> SetHueRange(0.66667, 0.0);
+  lut -> SetNumberOfColors(256);
+  lut -> SetRampToLinear();
+  lut -> Build();
+
+  
+	// map the contours to graphical primitives
+	vtkPolyDataMapper *contMapper = vtkPolyDataMapper::New();
+  contMapper->SetInput(normals -> GetOutput() );
+  contMapper->SetScalarRange(0,1);
+  contMapper->SetLookupTable(lut);
+  contMapper->ImmediateModeRenderingOn();
+
+ 
+ 	 // create an actor for the contours
+	vtkActor *contActor = vtkActor::New();
+	contActor->SetMapper(contMapper);
+ 
+	  // a renderer and render window
+	vtkRenderer *ren1 = vtkRenderer::New();
+	vtkRenderWindow *renWin = vtkRenderWindow::New();
+	renWin->AddRenderer(ren1);
+ 
+  	// an interactor
+	vtkRenderWindowInteractor *iren = vtkRenderWindowInteractor::New();
+	iren->SetRenderWindow(renWin);
+  
+  //The style
+  vtkInteractorStyleTerrain *terrain_style = vtkInteractorStyleTerrain::New();
+ 
+	// add the actors to the scene
+	ren1->SetBackground(0.6,0.6,0.6); 
+  ren1->GetActiveCamera()->SetViewUp(0,0,1);
+  ren1->GetActiveCamera()->SetParallelProjection(1);
+	ren1->AddActor(contActor);
+  iren -> SetInteractorStyle(terrain_style);
+  ren1->ResetCamera();
+
+ 
+  	// render an image (lights and cameras are created automatically)
+	renWin->Render();
+ 
+  	// begin mouse interaction
+	iren->Start();
+ 
+}
new file mode 100644
--- /dev/null
+++ b/vtk/terrain_test.tcl
@@ -0,0 +1,95 @@
+package require vtk
+package require vtkinteraction
+
+set Scale 5
+vtkLookupTable lut
+  lut SetHueRange 0.6 0
+  lut SetSaturationRange 1.0 0
+  lut SetValueRange 0.5 1.0
+
+set VTK_DATA_ROOT /usr/share/VTKData
+
+vtkDEMReader demModel
+  demModel SetFileName $VTK_DATA_ROOT/Data/SainteHelens.dem
+  demModel Update
+
+demModel Print
+
+set lo [expr $Scale * [lindex [demModel GetElevationBounds] 0]]
+set hi [expr $Scale * [lindex [demModel GetElevationBounds] 1]]
+
+vtkLODActor demActor
+
+# create a pipeline for each lod mapper
+vtkImageShrink3D shrink16
+  shrink16 SetShrinkFactors 16 16 1
+  shrink16 SetInputConnection [demModel GetOutputPort]
+  shrink16 AveragingOn
+
+vtkImageDataGeometryFilter geom16
+  geom16 SetInputConnection [shrink16 GetOutputPort]
+  geom16 ReleaseDataFlagOn
+
+vtkWarpScalar warp16
+  warp16 SetInputConnection [geom16 GetOutputPort]
+  warp16 SetNormal 0 0 1
+  warp16 UseNormalOn
+  warp16 SetScaleFactor $Scale
+  warp16 ReleaseDataFlagOn
+
+vtkElevationFilter elevation16
+  elevation16 SetInputConnection [warp16 GetOutputPort]
+  elevation16 SetLowPoint 0 0 $lo
+  elevation16 SetHighPoint 0 0 $hi
+  eval elevation16 SetScalarRange $lo $hi
+  elevation16 ReleaseDataFlagOn
+
+vtkPolyDataNormals normals16
+  normals16 SetInput [elevation16 GetPolyDataOutput]
+  normals16 SetFeatureAngle 60
+  normals16 ConsistencyOff
+  normals16 SplittingOff
+  normals16 ReleaseDataFlagOn
+
+vtkPolyDataMapper demMapper16
+  demMapper16 SetInputConnection [normals16 GetOutputPort]
+  eval demMapper16 SetScalarRange $lo $hi
+  demMapper16 SetLookupTable lut
+  demMapper16 ImmediateModeRenderingOn
+
+demMapper16 Update
+demActor AddLODMapper demMapper16
+
+# Create the RenderWindow, Renderer and both Actors
+#
+vtkRenderer ren1
+vtkRenderWindow renWin
+    renWin AddRenderer ren1
+vtkRenderWindowInteractor iren
+    iren SetRenderWindow renWin
+vtkInteractorStyleTerrain t
+iren SetInteractorStyle t
+ 
+# Add the actors to the renderer, set the background and size
+#
+ren1 AddActor demActor
+ren1 SetBackground .4 .4 .4
+
+iren AddObserver UserEvent {wm deiconify .vtkInteract}
+iren SetDesiredUpdateRate 1
+
+wm withdraw .
+proc TkCheckAbort {} {
+  set foo [renWin GetEventPending]
+    if {$foo != 0} {renWin SetAbortRender 1}
+}
+renWin AddObserver AbortCheckEvent {TkCheckAbort}
+
+[ren1 GetActiveCamera] SetViewUp 0 0 1
+[ren1 GetActiveCamera] SetPosition -99900 -21354 131801
+[ren1 GetActiveCamera] SetFocalPoint 41461 41461 2815
+ren1 ResetCamera
+[ren1 GetActiveCamera] Dolly 1.2
+ren1 ResetCameraClippingRange
+
+renWin Render