Fitting

Keywords fit and fitDiff provide the interpolated data and the derivative of the interpolated data respectively. This is useful when material laws, properties or some terms in variational formulation are given as a data file.

auto Xh = Pch<2>(mesh);
auto K = Xh->element();
K.on(_range=elements(mesh), _expr=fit(idv(T)[,dataFile(string),[type(string)]]));
Kd.on(_range=elements(mesh), _expr=fitDiff(idv(T)[,dataFile(string),[type(string)]]));

1. Options

--fit.datafile

the data file to interpolate (two columns)

--fit.type

P0, P1, Spline, Akima

--fit.P0

left = 0, right = 1, center = 2

--fit.P1_right

Kind of extention : zero = 0, constant = 1, extrapol = 2

--fit.P1_left

Kind of extention : zero = 0, constant = 1, extrapol = 2

--fit.Spline_right

natural = 0, clamped = 1

--fit.Spline_left

natural = 0, clamped = 1

    auto mesh = loadMesh(_mesh=new Mesh<Simplex<2>>);
    auto Xh = Pch<1>(mesh);
    auto T = Xh->element(); // the temperature (say)
    auto K = Xh->element(); // K(T) - the dependant of the temperature conductivity
    auto Kd= Xh->element(); // K'(T)
    auto T_K = Xh->element(); // T_ = true 
    auto T_Kd= Xh->element(); // 
    auto D_K = Xh->element(); // D_ = difference
    auto D_Kd= Xh->element(); // 
    T.on(_range=elements(mesh), _expr=Px());
    T_K.on(_range=elements(mesh),_expr=(5*Px()+sin(Px())));
    T_Kd.on(_range=elements(mesh),_expr=(5+cos(Px())));
    //double f(double t = 0.) { return 5. * t + sin(t); }
    auto f = [](double x = 0.) { return 5. * x + sin(x); };
#if 1
    auto e = exporter(_mesh = mesh );
#endif
    std::string datafilename = (fs::current_path()/"data.csv").string();
    if ( Environment::worldComm().isMasterRank() )
    {
        // Generates the datafile
        // we assume an unitsquare as mesh
        std::ofstream datafile( datafilename );
        datafile << "x, y\n";
        for(double t = -1; t < 2; t+=0.32)
            datafile << t << " , " << f(t) << "\n";
        datafile.close();
    }
    Environment::worldComm().barrier();

    std::vector<std::string> interpTypeRange = { "P0" , "P1", "Spline", "Akima" };
    for(int k = 0; k < interpTypeRange.size(); ++k )
    {
        std::string const& interpType = interpTypeRange[k];
        BOOST_TEST_MESSAGE( boost::format("test %1%")% interpType );
        // evaluate K(T) with the interpolation from the datafile
        K.on(_range=elements(mesh), _expr=fit( idv(T), datafilename, "x", "y", interpType ) );
        Kd.on(_range=elements(mesh), _expr=fitDiff( idv(T), datafilename, "x", "y", interpType ) );

        D_K.on(_range=elements(mesh),_expr=vf::abs(idv(K)-idv(T_K)));
        D_Kd.on(_range=elements(mesh),_expr=vf::abs(idv(Kd)-idv(T_Kd)));

        auto max_K = D_K.max();
        auto max_Kd= D_Kd.max();
#if 1
        e->step(k)->add("T",T);
        e->step(k)->add("K",K);
        e->step(k)->add("Kd",Kd);
        e->step(k)->add("T_K",T_K);
        e->step(k)->add("T_Kd",T_Kd);
        e->step(k)->add("D_K",D_K);
        e->step(k)->add("D_Kd",D_Kd);
        e->save();
#endif
        /// Note : the error has nothing to do with the mesh size but the step on the datafile
        switch( InterpolationTypeMap[interpType] )
        {
        case InterpolationType::P0: //P0 interpolation
        {
            BOOST_CHECK_SMALL(max_K, 0.95);
            BOOST_CHECK_SMALL(max_Kd, 6.0001); // the derivative is null
            break;
        }
        case InterpolationType::P1: // P1 interpolation
        {
            BOOST_CHECK_SMALL(max_K, 0.01);
            BOOST_CHECK_SMALL(max_Kd, 0.15);
            break;
        }
        case InterpolationType::Spline: // CSpline interpolation
        {
            BOOST_CHECK_SMALL(max_K, 0.01);
            BOOST_CHECK_SMALL(max_Kd, 0.15);
            break;
        }
        case InterpolationType::Akima: // Akima interpolation
        {
            BOOST_CHECK_SMALL(max_K, 0.016);
            BOOST_CHECK_SMALL(max_Kd, 0.03);
            break;
        }
        }

        BOOST_TEST_MESSAGE( boost::format("test %1% done")% interpType );
    }