adaptiveConvectionDiffusionReaction_example.cpp

Annotated source file

Here is the full file examples/adaptiveConvectionDiffusionReaction_example.cpp. Clicking on a function or class name will lead you to its reference documentation.

 
# include <gismo.h>
# include <gsAssembler/gsAdaptiveRefUtils.h>
 
using namespace std;
using namespace gismo;
 
int main(int argc, char *argv[])
{
   bool plot = false;
 
   gsCmdLine cmd("Example for solving a convection-diffusion problem.");
   cmd.addSwitch("plot", "Create a ParaView visualization file with the solution", plot);
   try { cmd.getValues(argc,argv); } catch (int rv) { return rv; }
 
   // --------------- specify exact solution and right-hand-side ---------------
 
   // Define exact solution (will be used for specifying Dirichlet boundary conditions
   //gsFunctionExpr<> g("if( y<-x/2-1/2, 1, 0 )", 2);
   gsFunctionExpr<> g("if( y>=0, if( x <=-1, 1, 0 ), 0 )", 2);
   // Define source function
   gsFunctionExpr<> rhs("0",2);
 
   // diffusion coefficient:
   gsFunctionExpr<> coeff_diff("0.000001","0","0","0.000001",2);
   // convection coefficient:
   gsFunctionExpr<> coeff_conv("3/sqrt(13)","-2/sqrt(13)",2);
   // reaction coefficient:
   gsFunctionExpr<> coeff_reac("0",2);
 
   // Print out source function and solution
   gsInfo<<"Source function " << rhs << "\n";
   gsInfo<<"Dirichlet boundary conditions "  << g << "\n\n";
 
 
   // --------------- read geometry from file ---------------
 
   // Read geometry from file
   // Read xml and create gsMultiPatch
   string fileSrc( "planar/lshape2d_3patches_thb.xml" );
   gsMultiPatch<real_t> patches;
   gsReadFile<real_t>( fileSrc, patches);
   gsInfo << "The domain is a "<< patches <<"\n";
 
   // Get all interfaces and boundaries:
   patches.computeTopology();
 
   // --------------- add bonudary conditions ---------------
   gsBoundaryConditions<> bcInfo;
 
   // For simplicity, set Dirichlet boundary conditions
   // given by exact solution g on all boundaries:
   for ( gsMultiPatch<>::const_biterator
            bit = patches.bBegin(); bit != patches.bEnd(); ++bit)
   {
       bcInfo.addCondition( *bit, condition_type::dirichlet, &g );
   }
 
   // --------------- define Pde ---------------
   gsConvDiffRePde<real_t> cdrPde(patches, bcInfo, & coeff_diff,& coeff_conv, & coeff_reac, & rhs);
 
 
   // --------------- set up basis ---------------
 
   // Copy basis from the geometry
   gsMultiBasis<> bases( patches );
 
 
   // Number of initial uniform refinement steps:
   int numInitUniformRefine  = 2;
 
   for (int i = 0; i < numInitUniformRefine; ++i)
     bases.uniformRefine();
 
 
   // --------------- set up adaptive refinement loop ---------------
 
   // Number of refinement loops to be done
   int numRefinementLoops = 3;
 
   // Specify cell-marking strategy...
   MarkingStrategy adaptRefCrit = PUCA;
   //MarkingStrategy adaptRefCrit = GARU;
   //MarkingStrategy adaptRefCrit = errorFraction;
 
   // ... and parameter.
   const real_t adaptRefParam = 0.7;
 
 
 
   // Construct assembler
   gsCDRAssembler<real_t> cdrAss( cdrPde, bases);
   // Set stabilization flag to 1 = SUPG
   cdrAss.options().setInt("Stabilization", stabilizerCDR::SUPG);
   // Compute Dirichlet values by L2-projection
   // Caution: Interpolation does not work for locally refined (T)HB-splines!
   cdrAss.options().setInt("DirichletValues",dirichlet::l2Projection);
 
   // --------------- adaptive refinement loop ---------------
 
   for( int refLoop = 0; refLoop <= numRefinementLoops; refLoop++)
   {
       gsInfo << "====== Loop " << refLoop << " of "
              <<numRefinementLoops<< " ======" << "\n";
 
       // --------------- solving ---------------
 
       // Generate system matrix and load vector
       cdrAss.assemble();
 
       // Solve the system
       gsMatrix<real_t> solVector =
           gsSparseSolver<>::BiCGSTABILUT( cdrAss.matrix() ).solve( cdrAss.rhs() );
 
       // Construct the solution as a scalar field
       gsField<> solField;
       solField = cdrAss.constructSolution(solVector);
 
       // --------------- error estimation/computation ---------------
 
       // Compute the H1-seminorm of the computed solution
       // ( which is, at least, equivalent to the energy norm in this example )
       // using the known exact solution.
       gsExprEvaluator<> ev;
       ev.setIntegrationElements(cdrAss.multiBasis());
       gsExprEvaluator<>::geometryMap Gm = ev.getMap(patches);
       gsExprEvaluator<>::variable is = ev.getVariable(solField.fields());
       auto ms = ev.getVariable(g, Gm);
 
       // Get the element-wise norms.
       ev.integralElWise( ( igrad(is,Gm) - igrad(ms)).sqNorm()*meas(Gm) );
       const std::vector<real_t> & eltErrs  = ev.elementwise();
 
       // --------------- adaptive refinement ---------------
 
       // Mark elements for refinement, based on the computed local errors and
       // the refinement-criterion and -parameter.
       std::vector<bool> elMarked( eltErrs.size() );
       gsMarkElementsForRef( eltErrs, adaptRefCrit, adaptRefParam, elMarked);
       gsInfo <<"Marked "<< std::count(elMarked.begin(), elMarked.end(), true) <<" elements.\n";
 
       // Refine the marked elements with a 1-ring of cells around marked elements
       gsRefineMarkedElements( cdrAss.multiBasis(), elMarked, 1 );
 
 
       // Call repair interfaces to make sure that the new meshes
       // match along patch interfaces.
       cdrAss.multiBasis().repairInterfaces( patches.interfaces() );
 
       cdrAss.refresh();
 
       // Export the final solution
       if( plot && refLoop == numRefinementLoops )
       {
           // Write the computed solution to paraview files
           gsWriteParaview<>( solField, "adaptRef", 1000, true);
       }
 
   }
 
   if( plot )
   {
       // Run paraview
       gsFileManager::open("adaptRef.pvd");
   }
   else
   {
       gsInfo<<"Done. No output created, re-run with --plot to get a ParaView "
               "file containing Plotting image data.\n";
   }
   return EXIT_SUCCESS;
 
}// end main