Work on non-homogeneous dirichlet bdry

examples/FEM4EllipticPDE_bhmag.cpp

         //Solver->SetSigmaFunction(implSigma);
         Solver->SetBoundaryStep(.05);
         Solver->SetGrid(MeshReader->GetOutput());
-        Solver->SetupPotential();
+        //Solver->SetupPotential();
         //Solver->SetGrid(rGrid);
         Solver->Solve(argv[3]);
 

examples/borehole/sphere.geo

  */
 
 radius = 2.25;     // Radius of the damn thing
-lc = radius/11;     //  0.25;   // Target element size
+lc = radius/3;     //  0.25;   // Target element size
 
 // Total Solution Space
 Box = 6*radius; // The down side of potential
 Z0 = -Box;
 Z1 =  Box;
 
-cellSize=radius/11; ///10;
+cellSize=radius/3; ///10;
 dd = 0 ; //  1e-5; //cellSize; // .01;
 pio2=Pi/2;
 

include/FEM4EllipticPDE.h

 #include "vtkCellData.h"
 
 #include <Eigen/IterativeLinearSolvers>
+#include <Eigen/SparseLU>
 #include "DCSurvey.h"
 
 namespace Lemma {

src/FEM4EllipticPDE.cpp

 
     void FEM4EllipticPDE::Solve( const std::string& resfile ) {
         ConstructAMatrix();
+        SetupPotential();
         //ConstructLoadVector();
 
         std::cout << "\nSolving" << std::endl;
+        std::cout << "   rows\tcols\n";
+        std::cout << "A: "  << A.rows() << "\t" << A.cols() << std::endl;
+        std::cout << "g: "  << g.rows() << "\t" << g.cols() << std::endl;
+
         ////////////////////////////////////////////////////////////
         // Solving:
         //Eigen::SimplicialCholesky<Eigen::SparseMatrix<Real>, Eigen::Lower > chol(A);  // performs a Cholesky factorization of A
         //VectorXr u = chol.solve(g);
 
+        //Eigen::SparseLU<Eigen::SparseMatrix<Real, Eigen::ColMajor>, Eigen::COLAMDOrdering<int> > solver;
+        //solver.analyzePattern(A);
+        //solver.factorize(A);
+        //VectorXr u = solver.solve(g);
+
+        //#ifdef CGSOLVE
         //Eigen::ConjugateGradient<Eigen::SparseMatrix<Real, Eigen::Lower >  Eigen::DiagonalPreconditioner > cg;
         Eigen::ConjugateGradient< Eigen::SparseMatrix<Real>  > cg(A);
-
         //Eigen::BiCGSTAB<Eigen::SparseMatrix<Real> > cg(A);
         //cg.setMaxIterations(3000);
         //cg.setTolerance(1e-28);
-
-        std::cout << "   rows\tcols\n";
-        std::cout << "A: "  << A.rows() << "\t" << A.cols() << std::endl;
-        std::cout << "g: "  << g.rows() << "\t" << g.cols() << std::endl;
-
         VectorXr u = cg.solve(g);
         std::cout << "#iterations:     " << cg.iterations() << std::endl;
         std::cout << "estimated error: " << cg.error()      << std::endl;
+        //#endif
 
         vtkDoubleArray *gArray = vtkDoubleArray::New();
         vtkDoubleArray *uArray = vtkDoubleArray::New();
             } else {
                 sigma_bar = 1.;
             }
+            sigma_bar = 1.;
 
-            auto W = VectorXr::Zero(4);
-            auto G = GradPhi.block<3,3>(1,1) * GradPhi.block<3,3>(1,1).transpose()*W;
 
             for (int ii=0; ii<4; ++ii) {
                 for (int jj=0; jj<4; ++jj) {
                     /* homogeneous Dirichlet boundary */
-                    if (jj == ii) {
+                    if (jj==ii && C(ii,3)==13.5) {
+                    //Real bb = vtkGrid->GetPointData()->GetScalars("HomogeneousDirichlet")->GetTuple(ID[ii])[0];
+                    //if (jj==ii && ((bb-1.)<1e-8) ) {
                         // Apply Homogeneous Dirichlet Boundary conditions
                         Real bb = vtkGrid->GetPointData()->GetScalars("HomogeneousDirichlet")->GetTuple(ID[ii])[0];
                         Real bdry = (1./(BndryH*BndryH))*BndrySigma*bb; // + sum;
                         //Real bdry = GradPhi.col(ii).tail<3>().dot(GradPhi.col(ii).tail<3>())*BndrySigma*bb; // + sum;
                         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 ) );
+                        //coeffs.push_back( Eigen::Triplet<Real> ( ID[ii], ID[jj], 1));
+                        //break;
                     } else {
                         coeffs.push_back( Eigen::Triplet<Real> ( ID[ii], ID[jj],        GradPhi.col(ii).tail<3>().dot(GradPhi.col(jj).tail<3>() ) * V * sigma_bar ) );
                     }
+                    /* */
+                    /* Inhomogeneous Dirichlet */
+                    //coeffs.push_back( Eigen::Triplet<Real> ( ID[ii], ID[jj], GradPhi.col(ii).tail<3>().dot(GradPhi.col(jj).tail<3>())*V*sigma_bar ));
                     // Stiffness matrix no longer contains boundary conditions...
                     //coeffs.push_back( Eigen::Triplet<Real> ( ID[ii], ID[jj], GradPhi.col(ii).tail<3>().dot(GradPhi.col(jj).tail<3>() ) * V * sigma_bar ) );
                 }
             }
-            std::cout <<  "\r" << (int)(1e2*((float)(ic) / (float)(vtkGrid->GetNumberOfCells()))) << std::flush ;
+            //std::cout <<  "\r" << (int)(1e2*((float)(ic) / (float)(vtkGrid->GetNumberOfCells()))) << std::flush ;
         }
         A.setFromTriplets(coeffs.begin(), coeffs.end());
         //std::cout << "A\n" << A << std::endl;
         std::cout << "\nBuilding load vector (g)" << std::endl;
         g = VectorXr::Zero(vtkGrid->GetNumberOfPoints());
         std::cout << "made g with "  << vtkGrid->GetNumberOfPoints() << " points" << std::endl;
-
+        VectorXr DB = VectorXr::Zero(vtkGrid->GetNumberOfPoints());
         for (int ic=0; ic < vtkGrid->GetNumberOfCells(); ++ic) {
 
             Eigen::Matrix<Real, 4, 4> C = Eigen::Matrix<Real, 4, 4>::Zero() ;
                 C(ii, 3) =  pts[2];
             }
 
-            Real V = (1./6.) * C.determinant();          // volume of tetrahedra
+            Real V = (1./6.) * C.determinant();              // volume of tetrahedra
+            Eigen::Matrix<Real, 4, 4> GradPhi = C.inverse(); // nabla \phi
 
             vtkIdList* Ids = vtkGrid->GetCell(ic)->GetPointIds();
             int ID[4];
                 ID[2] = Ids->GetId(2);
                 ID[3] = Ids->GetId(3);
 
-
+            /*
             Real avg(0);
             Real GG[4];
             for (int ii=0; ii<4; ++ii) {
             if ( std::abs( (GG[0]+GG[1]+GG[2]+GG[3])/4. - GG[0])  < 1e-5) {
                 avg = GG[0];
             }
+            */
+
+            VectorXr W = VectorXr::Zero(4);
+            for (int ii=0; ii<4; ++ii) {
+                W[ii] = vtkGrid->GetPointData()->GetScalars("HomogeneousDirichlet")->GetTuple(ID[ii])[0] *
+                        vtkGrid->GetPointData()->GetScalars("analytic_phi")->GetTuple(ID[ii])[0];
+                DB[ID[ii]] = vtkGrid->GetPointData()->GetScalars("HomogeneousDirichlet")->GetTuple(ID[ii])[0] *
+                             vtkGrid->GetPointData()->GetScalars("analytic_phi")->GetTuple(ID[ii])[0];
+            }
+            //auto G = GradPhi.block<3,4>(1,0).transpose()*GradPhi.block<3,4>(1,0)*W;
+            VectorXr G = GradPhi.block<3,4>(1,0).transpose()*GradPhi.block<3,4>(1,0)*W;
+            Real  sigma_bar(1.);
 
             for (int ii=0; ii<4; ++ii) {
             //     avg += vtkGrid->GetPointData()->GetScalars()->GetTuple(ID[ii])[0];
             //}
             //TODO this seems wrong!
             //avg /= 4.;
-                g(ID[ii]) +=  (V/4.) * ( vtkGrid->GetPointData()->GetScalars("G")->GetTuple(ID[ii])[0] )  ;
+                // TODO test code remove
+                g(ID[ii]) += V/4*(vtkGrid->GetPointData()->GetScalars("G")->GetTuple(ID[ii])[0])  ;
+                if (std::abs(C(ii,3)-13.5) > 1e-5) {
+                    g(ID[ii]) -= G[ii]*(V/3.)*sigma_bar;
+                }
                 //g(ID[ii]) +=  V/4*avg;
                   //g(ID[ii]) +=  6.67 *(V/4.) * avg;
             }
             //g(ID[0]) +=  (V/4.) * avg;
         }
+        //g -= A*DB;
 
     }
 
         //Eigen::SimplicialCholesky<Eigen::SparseMatrix<Real>, Eigen::Lower > chol(A);  // performs a Cholesky factorization of A
         //VectorXr u = chol.solve(g);
 
+        //Eigen::SparseLU<Eigen::SparseMatrix<Real, Eigen::ColMajor>, Eigen::COLAMDOrdering<int> > solver;
+        //solver.analyzePattern(A);
+        //solver.factorize(A);
+        //VectorXr u = solver.solve(g);
+
         //Eigen::ConjugateGradient<Eigen::SparseMatrix<Real, Eigen::Lower >  Eigen::DiagonalPreconditioner > cg;
         Eigen::ConjugateGradient< Eigen::SparseMatrix<Real>  > cg(A);
-
         //Eigen::BiCGSTAB<Eigen::SparseMatrix<Real> > cg(A);
         cg.setMaxIterations(3000);
         //cg.compute(A);
         //std::cout << "Computed     " << std::endl;
-
         VectorXr u = cg.solve(g);
         std::cout << "#iterations:     " << cg.iterations() << std::endl;
         std::cout << "estimated error: " << cg.error()      << std::endl;