Browse Source

Changed name of vtkValidPoint to HomogeneousDirichlet, as the name is more appropriate.

iss2
Trevor Irons 8 years ago
parent
commit
02db7f1d99
3 changed files with 9 additions and 41 deletions
  1. 1
    0
      examples/ResampleWithDataset.cpp
  2. 1
    34
      examples/borehole/sphere.sh
  3. 7
    7
      src/FEM4EllipticPDE.cpp

+ 1
- 0
examples/ResampleWithDataset.cpp View File

79
         Resample->SetInputData( uGrid );
79
         Resample->SetInputData( uGrid );
80
         Resample->SetSourceData( Stencil->GetOutput() );
80
         Resample->SetSourceData( Stencil->GetOutput() );
81
         Resample->SpatialMatchOn();
81
         Resample->SpatialMatchOn();
82
+        Resample->SetValidPointMaskArrayName("HomogeneousDirichlet");
82
         Resample->Update();
83
         Resample->Update();
83
 
84
 
84
     std::cout << *Resample << std::endl;
85
     std::cout << *Resample << std::endl;

+ 1
- 34
examples/borehole/sphere.sh View File

1
 #!/usr/bin/env bash
1
 #!/usr/bin/env bash
2
 gmsh  -3 -format vtk -o sphere.vtk  sphere.geo 
2
 gmsh  -3 -format vtk -o sphere.vtk  sphere.geo 
3
 gmsh  -2 -format stl -o sphereBox.stl  sphereBox.geo 
3
 gmsh  -2 -format stl -o sphereBox.stl  sphereBox.geo 
4
-//paraview --state=boundarySphere.pvsm 
5
 ../ResampleWithDataset  sphere.vtk  sphereBox.stl  MergedSphere.vtu 
4
 ../ResampleWithDataset  sphere.vtk  sphereBox.stl  MergedSphere.vtu 
6
-
7
-
8
 ../VTKEdgeGsphere MergedSphere.vtu   2.25  gravity  sphereGrav.vtu
5
 ../VTKEdgeGsphere MergedSphere.vtu   2.25  gravity  sphereGrav.vtu
9
 ../FEM4EllipticPDE_bhmag  sphereGrav.vtu  sphereGrav.vtu   sphereOutGrav.vtu 
6
 ../FEM4EllipticPDE_bhmag  sphereGrav.vtu  sphereGrav.vtu   sphereOutGrav.vtu 
10
 
7
 
11
-paraview --state=sliceGravSphere3.pvsm
12
-
13
-#IMPORTANT! 
14
-#   Select ResampleWithDataset1
15
-#	-->File -> SaveData -> MergedSphere.vtu 
16
-
17
-# Set up problem 
18
-
19
-
20
-# Gravity problem
21
-#../VTKEdgeGsphere MergedSphere.vtu   2.25  gravity  sphereGrav.vtu
22
-#../FEM4EllipticPDE_bhmag  sphereGrav.vtu  sphereGrav.vtu   sphereOutGrav.vtu 
23
-
24
-#paraview --state=sliceSphere.pvsm
25
-
26
-
27
-# --data=sphereOut.vtu
28
-
29
-# even finer meshing
30
-# radius = 4, u = \pm .9              analytic = \pm 1.33
31
-
32
-# Finer meshing
33
-# radius = 5, u = \pm 1.65            analytic = \pm 1.67
34
-# radius = 4, u = \pm 1.09            analytic = \pm 1.33
35
-
36
-# radius = 5, u = \pm 2.5             analytic = \pm 1.67
37
-# radius = 4, u = \pm 1.66            analytic = \pm 1.33
38
-# radius = 3, u= \pm .926 -.985       analytic = \pm 1
39
-# radius = 2, u = \pm.42              analytic = \pm .667
40
-# for radius=1,  u = \pm .111         analytic = \pm .333 
41
-# for radius=.5, u = \pm .0477 -.0462 analytic = \pm .167
8
+#paraview --data=sphereOutGrav.vtu 

+ 7
- 7
src/FEM4EllipticPDE.cpp View File

199
 
199
 
200
 
200
 
201
 
201
 
202
-    void FEM4EllipticPDE::Solve( const std::string& resfile) {
202
+    void FEM4EllipticPDE::Solve( const std::string& resfile ) {
203
         ConstructAMatrix();
203
         ConstructAMatrix();
204
         //ConstructLoadVector();
204
         //ConstructLoadVector();
205
 
205
 
213
         Eigen::ConjugateGradient< Eigen::SparseMatrix<Real>  > cg(A);
213
         Eigen::ConjugateGradient< Eigen::SparseMatrix<Real>  > cg(A);
214
 
214
 
215
         //Eigen::BiCGSTAB<Eigen::SparseMatrix<Real> > cg(A);
215
         //Eigen::BiCGSTAB<Eigen::SparseMatrix<Real> > cg(A);
216
-        cg.setMaxIterations(3000);
217
-        cg.setTolerance(1e-28);
216
+        //cg.setMaxIterations(3000);
217
+        //cg.setTolerance(1e-28);
218
 
218
 
219
-        std::cout << "  rows\tcols\n";
219
+        std::cout << "   rows\tcols\n";
220
         std::cout << "A: "  << A.rows() << "\t" << A.cols() << std::endl;
220
         std::cout << "A: "  << A.rows() << "\t" << A.cols() << std::endl;
221
         std::cout << "g: "  << g.rows() << "\t" << g.cols() << std::endl;
221
         std::cout << "g: "  << g.rows() << "\t" << g.cols() << std::endl;
222
 
222
 
265
         A.resize(vtkGrid->GetNumberOfPoints(), vtkGrid->GetNumberOfPoints());
265
         A.resize(vtkGrid->GetNumberOfPoints(), vtkGrid->GetNumberOfPoints());
266
         std::vector< Eigen::Triplet<Real> > coeffs;
266
         std::vector< Eigen::Triplet<Real> > coeffs;
267
 
267
 
268
-        if ( !vtkGrid->GetPointData()->GetScalars("vtkValidPointMask") ) {
269
-            throw std::runtime_error("No vtkValidPointMask");
268
+        if ( !vtkGrid->GetPointData()->GetScalars("HomogeneousDirichlet") ) {
269
+            throw std::runtime_error("No HomogeneousDirichlet boundary conditions in input file.");
270
         }
270
         }
271
 
271
 
272
         if ( !vtkGrid->GetCellData()->GetScalars("G") && !vtkGrid->GetPointData()->GetScalars("G") ) {
272
         if ( !vtkGrid->GetCellData()->GetScalars("G") && !vtkGrid->GetPointData()->GetScalars("G") ) {
325
                     /* homogeneous Dirichlet boundary */
325
                     /* homogeneous Dirichlet boundary */
326
                     if (jj == ii) {
326
                     if (jj == ii) {
327
                         // Apply Homogeneous Dirichlet Boundary conditions
327
                         // Apply Homogeneous Dirichlet Boundary conditions
328
-                        Real bb = vtkGrid->GetPointData()->GetScalars("vtkValidPointMask")->GetTuple(ID[ii])[0];
328
+                        Real bb = vtkGrid->GetPointData()->GetScalars("HomogeneousDirichlet")->GetTuple(ID[ii])[0];
329
                         Real bdry = (1./(BndryH*BndryH))*BndrySigma*bb; // + sum;
329
                         Real bdry = (1./(BndryH*BndryH))*BndrySigma*bb; // + sum;
330
                         //Real bdry = GradPhi.col(ii).tail<3>().dot(GradPhi.col(ii).tail<3>())*BndrySigma*bb; // + sum;
330
                         //Real bdry = GradPhi.col(ii).tail<3>().dot(GradPhi.col(ii).tail<3>())*BndrySigma*bb; // + sum;
331
                         coeffs.push_back( Eigen::Triplet<Real> ( ID[ii], ID[jj], bdry + GradPhi.col(ii).tail<3>().dot(GradPhi.col(jj).tail<3>() ) * V * sigma_bar ) );
331
                         coeffs.push_back( Eigen::Triplet<Real> ( ID[ii], ID[jj], bdry + GradPhi.col(ii).tail<3>().dot(GradPhi.col(jj).tail<3>() ) * V * sigma_bar ) );

Loading…
Cancel
Save