a_fcriteria_coulomb.cxx

Go to the documentation of this file.

/*
Copyright 2016 Pierre SMARS (smars@yuntech.edu.tw)
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 2 of the License, or
(at your option) any later version.
 
This program 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 General Public License for more details.
 
You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/
#include "a_fcriteria_coulomb.h"
#include "a_face_2d4.h"
 
#include "config.h"
 
//---------------------------------------------------------------------------
const std::string a_fcriteria_coulomb::help()
{
        std::ostringstream o;
        o << a_fcriteria::help();
        o << "********************" << std::endl;
        o << "a_fcriteria_coulomb:" << std::endl;
        o << "********************" << std::endl;
        o << "This class derives from a _fcriteria_in," << std::endl;
        o << "it implements a specific resistance criterion:" << std::endl;
        o << "The excentricity of the resultant of the force on a face is required to remain in between specific bounds" << std::endl;
        o << "And the shear force on a face needs to respect Coulomb criteria: |Q| <= f |N| + c" << std::endl;
        o << "With |Q| = |Qx + Qy|" << std::endl;
        o << "f 'value': set the bound (default: 0.5 = 26.6 deg.)" << std::endl;
        o << "f: get the bound" << std::endl;
        o << "c 'value': set the cohesion (default: 0.)" << std::endl;
        o << "c: get the cohesion" << std::endl;
        return o.str();
}
//---------------------------------------------------------------------------
bool a_fcriteria_coulomb::ok(const a_face * f0) const
{
        if (!a_fcriteria_in::ok(f0))
                return false;
        const a_face_2d4 * f = cf_2d4(f0);
        double N = f->N();
        if (N>0.)
                return false;
        double Qx = f->Qx();
        double Qy = f->Qy();
        double Q = sqrt(Qx*Qx+Qy*Qy);
        if (Q>-f_*N+c_)
                return false;
        return true;
}
//---------------------------------------------------------------------------
double a_fcriteria_coulomb::penalty(const a_face * f0) const
{
        double p1 = a_fcriteria_in::penalty(f0);
        const a_face_2d4 * f = cf_2d4(f0);
        double N = f->N();
        double Qx = f->Qx();
        double Qy = f->Qy();
        double Q = sqrt(Qx*Qx+Qy*Qy);
        if (Q<=f_*N+c_)
                return p1;
        double p2 = Q-f_*N-c_;
        if (N!=0)
        {
                p2 /= N;
                return sqrt(p1*p1+p2*p2);
        }
        return c_/Q;
}
//---------------------------------------------------------------------------