Homec4science
Diffusion tamaas (maxwell_viscoelastic-branch) Browse

tamaas/src/function.cpp
7d473b01cbbamaxwell_viscoelastic…

Last Change
Raw File

function.cpp
View Options

/**
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @section LICENSE
*
* Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de
* Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des
* Solides)
*
* Tamaas is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Tamaas is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
#include "function.hh"
/* -------------------------------------------------------------------------- */
INTER_FUNC::INTER_FUNC() {}
/* ------------------------------------------------------------------------ */
INTER_FUNC::~INTER_FUNC() {}
/* ------------------------------------------------------------------------ */
INTER_FUNC::INTER_FUNC(INTER_FUNCTION_TYPE type_, std::vector<Real> points_, Real a_)
{
type = type_;
sm_error = false;
points.resize(points_.size());
for (int i = 0; i < int(points.size()); i++)
points[i] = points_[i];
if (type == SQ_ROOT)
{
if (points.size() != 2)
{
std::cout << std::endl << "ERROR 1a: something wrong with interpolation; coef[0] = " << coef[0] << std::endl;
exit(1);
}
coef.resize(3);
coef[2] = a_;
Real y1_sq = points[0]*points[0],
y2_sq = points[1]*points[1],
c0 = y2_sq/(y1_sq-y2_sq),
c1 = (2+c0),
c2 = (1 - coef[2]*coef[2] +0.75*c0),
x01 = 0.5*(c1+sqrt(c1*c1-4*c2)),
x02 = 0.5*(c1-sqrt(c1*c1-4*c2));
if (y1_sq <= y2_sq)
{
if (x01 > 0.75 && x02 <= 0.75)
{
coef[0] = sqrt((y1_sq-y2_sq)/(0.75-x01));
coef[1] = x01;
}
else if (x02 > 0.75 && x01 <= 0.75)
{
coef[0] = sqrt((y1_sq-y2_sq)/(0.75-x02));
coef[1] = x02;
}
else if (x01 > 0.75 && x02 <= 0.75 && a_ <= 0.5)
{
coef[0] = sqrt((y1_sq-y2_sq)/(0.75-x02));
coef[1] = x02;
}
else if (x02 > 0.75 && x01 <= 0.75 && a_ <= 0.5)
{
coef[0] = sqrt((y1_sq-y2_sq)/(0.75-x01));
coef[1] = x01;
}
else
{
std::cout << std::endl << "ERROR 1b: something wrong with interpolation; coef[0] = " << coef[0] << std::endl;
exit(1);
}
}
else
{
coef[0] = 0;
coef[1] = 0;
coef[2] = 0;
sm_error = true;
}
return;
}
else
{
std::cout << std::endl << "ERROR 3: wrong constructor" << std::endl;
exit(1);
}
}
/* ------------------------------------------------------------------------ */
INTER_FUNC::INTER_FUNC(INTER_FUNCTION_TYPE type_, std::vector<Real> points_)
{
type = type_;
sm_error = false;
points.resize(points_.size());
for (int i = 0; i < int(points.size()); i++)
points[i] = points_[i];
if (type == LINEAR)
{
if (points.size() != 2)
{
std::cout << "ERROR 21: Wrong initialization" << std::endl;
exit(1);
}
coef.resize(2);
coef[0] = points[0];
coef[1] = points[1]-points[0];
return;
}
else if (type == PARABOLA)
{
if (points.size() != 3)
{
std::cout << "ERROR 22: Wrong initialization" << std::endl;
exit(1);
}
coef.resize(3);
coef[2] = points[0];
coef[1] = 4*points[1]-3*points[0]-points[2];
coef[0] = 4*points[1]-4*points[0]-2*coef[1];
return;
}
else if (type != ZERO)
{
std::cout << std::endl << "ERROR 4: wrong constructor" << std::endl;
exit(1);
}
}
/* ------------------------------------------------------------------------ */
Real INTER_FUNC::get_value(Real x)
{
if (type == ZERO)
{
return 0.;
}
else if (type == LINEAR)
{
return coef[0] + coef[1]*x;
}
else if (type == SQ_ROOT)
{
Real a_sq = coef[2]*coef[2],
x_x0_sq = (x-coef[1])*(x-coef[1]);
if (a_sq > x_x0_sq)
return coef[0]*sqrt(a_sq - x_x0_sq);
else
return 0.;
}
else if (type == PARABOLA)
{
return coef[0]*x*x + coef[1]*x + coef[2];
}
std::cout << std::endl << "ERROR 4.5: unknown interpolation type" << std::endl;
exit(1);
return 0.;
}
/* ------------------------------------------------------------------------ */
Real INTER_FUNC::get_derivative(Real x)
{
if (type == ZERO)
{
return 0.;
}
else if (type == LINEAR)
{
return coef[1]*x;
}
else if (type == SQ_ROOT)
{
Real a_sq = coef[2]*coef[2],
x_x0_sq = (x-coef[1])*(x-coef[1]);
if (a_sq > x_x0_sq)
return coef[0]*(coef[1]-x)/sqrt(a_sq - x_x0_sq);
else
return 0.;
}
else if (type == PARABOLA)
{
return 2*coef[0]*x + coef[1];
}
std::cout << std::endl << "ERROR 3: something wrong with interpolation" << std::endl;
exit(1);
}
c4science · Help