#include <feel/feelmodels/maxwell/maxwell.hpp>
#include <feel/feelmodels/modelcore/options.hpp>
#include <feel/feelfilters/loadmesh.hpp>
using namespace Feel;
void printCvg( std::ostream& out, std::vector<double> h, std::vector<std::vector<double>> err )
{
auto errA = err[0];
auto errRelA = err[1];
auto errB = err[2];
auto errRelB = err[3];
boost::format fmter("%1% %|14t|");
boost::format fmter_endl("%1%\n");
out << fmter % "h";
out << fmter % "errA_L2";
out << fmter % "errA_Curl";
out << fmter % "errRelA_L2";
out << fmter % "errRelA_Curl";
out << fmter % "errB_L2";
out << fmter % "errB_Div";
out << fmter % "errRelB_L2";
out << fmter_endl % "errRelB_Div";
int maxiter = h.size();
for( int i = 0; i < maxiter; ++i )
{
out << fmter % h[i];
out << fmter % errA[2*i];
out << fmter % errA[2*i+1];
out << fmter % errRelA[2*i];
out << fmter % errRelA[2*i+1];
out << fmter % errB[2*i];
out << fmter % errB[2*i+1];
out << fmter % errRelB[2*i];
out << fmter_endl % errRelB[2*i+1];
}
}
// template <int OrderT>
void
runApplicationMaxwell()
{
using model_type = FeelModels::Maxwell< Simplex<FEELPP_DIM> >;
using mesh_type = typename model_type::mesh_type;
using exporter_type = Exporter<mesh_type>;
using exporter_ptrtype = std::shared_ptr<exporter_type>;
std::shared_ptr<model_type> maxwell( new model_type("maxwell", false) );
int maxiter = ioption("cvg.nb-iter");
double factor = doption("cvg.factor");
double size = doption("gmsh.hsize");
exporter_ptrtype e( exporter_type::New("convergence") );
std::vector<double> h(maxiter);
std::vector<double> normA(2*maxiter), errA(2*maxiter), errRelA(2*maxiter);
std::vector<double> normB(2*maxiter), errB(2*maxiter), errRelB(2*maxiter);
for( int i = 0; i < maxiter; ++i )
{
h[i] = size;
auto mesh = loadMesh( _mesh=new mesh_type, _h=size );
maxwell->setMesh(mesh);
maxwell->init();
maxwell->printAndSaveInfo();
maxwell->solve();
auto functions = maxwell->modelProperties().functions();
auto params = maxwell->modelProperties().parameters().toParameterValues();
functions.setParameterValues( params );
e->step(i)->setMesh(mesh);
auto A_h = maxwell->fieldMagneticPotential();
auto AExpr = expr<FEELPP_DIM,1>(functions["A_exact"].expressionString() );
AExpr.setParameterValues( params );
auto Ah = maxwell->spaceMagneticPotential();
auto A_ex = Ah->element(AExpr);
normA[2*i] = normL2( _range=Ah->template rangeElements<0>(),
_expr=idv(A_h) );
errA[2*i] = normL2( _range=Ah->template rangeElements<0>(),
_expr=idv(A_h)-idv(A_ex) );
errRelA[2*i] = errA[2*i]/normA[2*i];
normA[2*i+1] = normL2( _range=Ah->template rangeElements<0>(),
_expr=idv(A_h) );
normA[2*i+1] += normL2( _range=Ah->template rangeElements<0>(),
_expr=curlv(A_h) );
errA[2*i+1] = normL2( _range=Ah->template rangeElements<0>(),
_expr=idv(A_h)-idv(A_ex) );
errA[2*i+1] += normL2( _range=Ah->template rangeElements<0>(),
_expr=curlv(A_h)-curlv(A_ex) );
errRelA[2*i+1] = errA[2*i+1]/normA[2*i+1];
e->step(i)->add("A_h", "A_h", A_h);
e->step(i)->add("A_ex", "A_ex", A_ex);
auto B_h = maxwell->fieldMagneticField();
#if FEELPP_DIM==3
auto BExpr = expr<FEELPP_DIM,1>(functions["B_exact"].expressionString() );
#else
auto BExpr = expr(functions["B_exact"].expressionString() );
#endif
BExpr.setParameterValues( params );
auto Bh = maxwell->spaceMagneticField();
auto B_ex = Bh->element(BExpr);
normB[2*i] = normL2( _range=Bh->template rangeElements<0>(),
_expr=idv(B_h) );
errB[2*i] = normL2( _range=Bh->template rangeElements<0>(),
_expr=idv(B_h)-idv(B_ex) );
errRelB[2*i] = errB[2*i]/normB[2*i];
#if FEELPP_DIM==3
normB[2*i+1] = normL2( _range=Bh->template rangeElements<0>(),
_expr=idv(B_h) );
normB[2*i+1] += normL2( _range=Bh->template rangeElements<0>(),
_expr=divv(B_h) );
errB[2*i+1] = normL2( _range=Bh->template rangeElements<0>(),
_expr=idv(B_h)-idv(B_ex) );
errB[2*i+1] += normL2( _range=Bh->template rangeElements<0>(),
_expr=divv(B_h)-divv(B_ex) );
errRelB[2*i+1] = errB[2*i+1]/normB[2*i+1];
#endif
e->step(i)->add("B_h", "B_h", B_h);
e->step(i)->add("B_ex", "N_ex", B_ex);
e->save();
size /= factor;
}
std::vector<std::vector<double>> err = {errA, errRelA, errB, errRelB};
Feel::cout << std::endl;
if( Environment::isMasterRank() )
{
printCvg( std::cout, h, err );
std::ofstream file ( soption("cvg.filename") );
if( file )
{
printCvg( file, h, err );
file.close();
}
}
}
int main(int argc, char** argv) {
po::options_description maxwelloptions( "application maxwell options" );
maxwelloptions.add( toolboxes_options("maxwell") );
maxwelloptions.add_options()
("cvg.nb-iter", po::value<int>()->default_value(1), "number of convergence iteration")
("cvg.factor", po::value<double>()->default_value(2), "factor of mesh size by iteration")
("cvg.filename", po::value<std::string>()->default_value("cvg.dat"), "filename to store convergence results")
;
Environment env( _argc=argc, _argv=argv,
_desc=maxwelloptions,
_about=about(_name="cvg_maxwell",
_author="Feel++ Consortium",
_email="feelpp-devel@feelpp.org"));
runApplicationMaxwell();
return 0;
}