Adding Gilbert Strang approach to Dirichlet conditions.

examples/borehole/sphere.geo

@@ -8,10 +8,10 @@
  */
 
 radius = 2.25;     // Radius of the damn thing
-lc = radius/3;     //  0.25;   // Target element size
+lc = radius/4;     //  0.25;   // Target element size
 
 // Total Solution Space
-Box = 6*radius; // The down side of potential
+Box = 3*radius; // The down side of potential
 X0 = -Box;
 X1 =  Box;
 Y0 = -Box;
@@ -19,7 +19,7 @@ Y1 =  Box;
 Z0 = -Box;
 Z1 =  Box;
 
-cellSize=radius/3; ///10;
+cellSize=radius/32; ///10;
 dd = 0 ; //  1e-5; //cellSize; // .01;
 pio2=Pi/2;
 

examples/borehole/sphereBox.geo

@@ -11,10 +11,10 @@ D0 = 10;          // Top of magnet
 D1 = 11;          // Bottom of magnet
 radius = 2.25;     // Radius of the damn thing
 
-lc = 2*radius;   //  0.25;   // Target element size
+lc = radius;   //  0.25;   // Target element size
 
 // Total Solution Space
-Box = 6*radius; // The down side of potential
+Box = 3*radius; // The down side of potential
 X0 = -Box;
 X1 =  Box;
 Y0 = -Box;

src/FEM4EllipticPDE.cpp

@@ -221,8 +221,8 @@ namespace Lemma {
 
         //#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);
+        //Eigen::ConjugateGradient< Eigen::SparseMatrix<Real>  > cg(A);
+        Eigen::BiCGSTAB<Eigen::SparseMatrix<Real> > cg(A);
         //cg.setMaxIterations(3000);
         //cg.setTolerance(1e-28);
         VectorXr u = cg.solve(g);
@@ -340,21 +340,24 @@ namespace Lemma {
 
 
             for (int ii=0; ii<4; ++ii) {
-                for (int jj=0; jj<4; ++jj) {
-                    /* homogeneous Dirichlet boundary */
-                    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 {
+                int bbi = vtkGrid->GetPointData()->GetScalars("HomogeneousDirichlet")->GetTuple(ID[ii])[0];
+                if (bbi) {
+                    coeffs.push_back( Eigen::Triplet<Real> ( ID[ii], ID[ii], 1));
+                } else {
+                    for (int jj=0; jj<4; ++jj) {
+                        /* homogeneous Dirichlet boundary */
+                        //if (jj==ii && C(ii,3)==13.5) {
+                        //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 ));
@@ -362,6 +365,7 @@ namespace Lemma {
                     //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 ;
         }
         A.setFromTriplets(coeffs.begin(), coeffs.end());
@@ -409,6 +413,7 @@ namespace Lemma {
             }
             */
 
+            /*
             VectorXr W = VectorXr::Zero(4);
             for (int ii=0; ii<4; ++ii) {
                 W[ii] = vtkGrid->GetPointData()->GetScalars("HomogeneousDirichlet")->GetTuple(ID[ii])[0] *
@@ -419,6 +424,7 @@ namespace Lemma {
             //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];
@@ -427,10 +433,14 @@ namespace Lemma {
             //TODO this seems wrong!
             //avg /= 4.;
                 // 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;
+                if (vtkGrid->GetPointData()->GetScalars("HomogeneousDirichlet")->GetTuple(ID[ii])[0]) {
+                    g(ID[ii]) += vtkGrid->GetPointData()->GetScalars("analytic_phi")->GetTuple(ID[ii])[0];
+                } else {
+                    g(ID[ii]) += (3.0*V)*(vtkGrid->GetPointData()->GetScalars("G")->GetTuple(ID[ii])[0]); // Why 3.0??
                 }
+                //if (std::abs(C(ii,3)-13.5) > 1e-5) {
+                //    g(ID[ii]) -= G[ii]*(V)*sigma_bar;
+                //}
                 //g(ID[ii]) +=  V/4*avg;
                   //g(ID[ii]) +=  6.67 *(V/4.) * avg;
             }