Skip to content
Snippets Groups Projects
Select Git revision
  • dfb7dee1a8f2a3392ff59829c20a30c62096f7c4
  • master default protected
  • alphashapes
  • quadMeshingTools
  • cygwin_conv_path
  • macos_arm64
  • add-transfiniteautomatic-to-geo
  • patch_releases_4_10
  • HierarchicalHDiv
  • isuruf-master-patch-63355
  • hyperbolic
  • hexdom
  • hxt_update
  • jf
  • 1618-pythonocc-and-gmsh-api-integration
  • octreeSizeField
  • hexbl
  • alignIrregularVertices
  • getEdges
  • patch_releases_4_8
  • isuruf-master-patch-51992
  • gmsh_4_11_0
  • gmsh_4_10_5
  • gmsh_4_10_4
  • gmsh_4_10_3
  • gmsh_4_10_2
  • gmsh_4_10_1
  • gmsh_4_10_0
  • gmsh_4_9_5
  • gmsh_4_9_4
  • gmsh_4_9_3
  • gmsh_4_9_2
  • gmsh_4_9_1
  • gmsh_4_9_0
  • gmsh_4_8_4
  • gmsh_4_8_3
  • gmsh_4_8_2
  • gmsh_4_8_1
  • gmsh_4_8_0
  • gmsh_4_7_1
  • gmsh_4_7_0
41 results

meshGFaceOptimize.cpp

Blame
  • Forked from gmsh / gmsh
    Source project has a limited visibility.
    Polynomial.cpp 13.63 KiB
    #include <cmath>
    #include <sstream>
    #include <stack>
    
    #include "Polynomial.h"
    
    using namespace std;
    
    const char Polynomial::coefName[3] = {'x', 'y', 'z'};
    
    Polynomial::Polynomial(double coef,
                           int powerX,
                           int powerY,
                           int powerZ){
      nMon = 1;
      mon  = new monomial_t[1];
    
      mon[0].coef     = coef;
      mon[0].power[0] = powerX;
      mon[0].power[1] = powerY;
      mon[0].power[2] = powerZ;
    }
    
    Polynomial::Polynomial(const Polynomial& other){
      nMon = other.nMon;
      mon  = copyMonomial(other.mon, nMon);
    }
    
    Polynomial::Polynomial(void){
      nMon = 0;
      mon  = NULL;
    }
    
    Polynomial::~Polynomial(void){
      if(mon)
        delete[] mon;
    }
    
    void Polynomial::derivative(int dim){
      // Take derivative //
      for(int i = 0; i < nMon; i++){
        mon[i].coef *= mon[i].power[dim];
        mon[i].power[dim] -= 1;
      }
    
      // Remove zero monomials //
      int N = 0;
      stack<monomial_t*> s;
    
      for(int i = 0; i < nMon; i++){
        if(mon[i].coef != 0.0){
          s.push(&mon[i]);
          N++;
        }
      }
    
      // If no monomial any more ---> return zero polynomial
      if(!N){
        delete[] mon;
    
        mon  = zeroPolynomial();
        nMon = 1;
        return;
      }
    
      // If no zero found ---> return;
      if(N == nMon)
        return;
    
      // Else, remove them //
      monomial_t* tmp = mon;
    
      mon  = new monomial_t[N];
      nMon = N;
    
      for(int i = N - 1; i >= 0; i--){
        mon[i] = *(s.top());
        s.pop();
      }
    
      delete[] tmp;
    
      // Sort resulting monomial and return //
      sort(mon, nMon);
    
      return;
    }
    
    vector<Polynomial> Polynomial::gradient(void) const{
      vector<Polynomial> grad(3);
    
      // Copy Polynomial //
      grad[0] = *this;
      grad[1] = *this;
      grad[2] = *this;
    
      // Derivative with respect to each direction //
      grad[0].derivative(0);
      grad[1].derivative(1);
      grad[2].derivative(2);
    
      return grad;
    }
    
    vector<Polynomial> Polynomial::curl(const vector<Polynomial>& p){
      vector<Polynomial> rot(3);
    
      // Partial Derivatives //
      Polynomial dP0d1 = p[0];
      Polynomial dP0d2 = p[0];
      Polynomial dP1d0 = p[1];
      Polynomial dP1d2 = p[1];
      Polynomial dP2d0 = p[2];
      Polynomial dP2d1 = p[2];
    
      dP0d1.derivative(1);
      dP0d2.derivative(2);
      dP1d0.derivative(0);
      dP1d2.derivative(2);
      dP2d0.derivative(0);
      dP2d1.derivative(1);
    
      // Curl //
      rot[0] = dP2d1 - dP1d2;
      rot[1] = dP0d2 - dP2d0;
      rot[2] = dP1d0 - dP0d1;
    
      // Return //
      return rot;
    }
    
    Polynomial Polynomial::divergence(const vector<Polynomial>& p){
      // Partial Derivatives //
      Polynomial dP0d0 = p[0];
      Polynomial dP1d1 = p[1];
      Polynomial dP2d2 = p[2];
    
      dP0d0.derivative(0);
      dP1d1.derivative(1);
      dP2d2.derivative(2);
    
      // Return Div //
      return dP0d0 + dP1d1 + dP2d2;
    }
    
    double Polynomial::at
      (double x, double y, double z) const{
    
      double val = 0;
      for(int i = 0; i < nMon; i++){
        val += mon[i].coef * pow(x, mon[i].power[0])
                           * pow(y, mon[i].power[1])
                           * pow(z, mon[i].power[2]);
      }
    
      return val;
    }
    
    fullVector<double> Polynomial::at(const std::vector<Polynomial>& P,
                                      double x,
                                      double y,
                                      double z){
      fullVector<double> val(3);
    
      val(0) = P[0].at(x, y, z);
      val(1) = P[1].at(x, y, z);
      val(2) = P[2].at(x, y, z);
    
      return val;
    }
    
    
    Polynomial Polynomial::operator+(const Polynomial& other) const{
      Polynomial newP;
    
      newP.nMon = mergeMon(mon, nMon, other.mon, other.nMon, &newP.mon);
    
      return newP;
    }
    
    Polynomial Polynomial::operator-(const Polynomial& other) const{
      const int   otherNMon  = other.nMon;
      monomial_t* otherMinus = copyMonomial(other.mon, otherNMon);
    
      Polynomial  newP;
    
      mult(otherMinus, otherNMon, -1);
      newP.nMon = mergeMon(mon, nMon, otherMinus, otherNMon, &newP.mon);
    
      delete[] otherMinus;
      return newP;
    }
    
    Polynomial Polynomial::operator*(const Polynomial& other) const{
      Polynomial newP;
    
      newP.nMon = mult(mon, nMon, other.mon, other.nMon, &newP.mon);
    
      return newP;
    }
    
    Polynomial Polynomial::operator*(double alpha) const{
      Polynomial newP;
    
      newP.mon  = copyMonomial(mon, nMon);
      newP.nMon = nMon;
    
      newP.mul(alpha);
    
      return newP;
    }
    
    
    void Polynomial::add(const Polynomial& other){
      monomial_t* tmp = mon;
    
      nMon = mergeMon(mon, nMon, other.mon, other.nMon, &mon);
    
      delete[] tmp;
    }
    
    void Polynomial::sub(const Polynomial& other){
      const int   otherNMon  = other.nMon;
      monomial_t* otherMinus = copyMonomial(other.mon, otherNMon);
      monomial_t* tmp        = mon;
    
      mult(otherMinus, otherNMon, -1);
      nMon = mergeMon(mon, nMon, otherMinus, otherNMon, &mon);
    
      delete[] otherMinus;
      delete[] tmp;
    }
    
    void Polynomial::mul(const Polynomial& other){
      monomial_t* tmp = mon;
    
      nMon = mult(mon, nMon, other.mon, other.nMon, &mon);
    
      delete[] tmp;
    }
    
    void Polynomial::mul(double alpha){
      for(int i = 0; i < nMon; i++)
        mon[i].coef *= alpha;
    }
    
    void Polynomial::power(int n){
      if (n < 0)
        return;
    
      switch(n){
      case 0:
        delete[] mon;
    
        mon  = unitPolynomial();
        nMon = 1;
    
        break;
    
      case 1:
        break;
    
      default:
        Polynomial old = *this;
    
        for(int i = 1; i < n; i++)
          mul(old);
    
        break;
      }
    }
    
    Polynomial Polynomial::compose(const Polynomial& other) const{
      stack<monomial_t> stk;
    
      for(int i = 0; i < nMon; i++)
        compose(&mon[i], other, &stk);
    
      return polynomialFromStack(stk);
    }
    
    Polynomial Polynomial::compose(const Polynomial& otherA,
                                   const Polynomial& otherB) const{
      stack<monomial_t> stk;
    
      for(int i = 0; i < nMon; i++)
        compose(&mon[i], otherA, otherB, &stk);
    
      return polynomialFromStack(stk);
    }
    
    Polynomial Polynomial::compose(const Polynomial& otherA,
                                   const Polynomial& otherB,
                                   const Polynomial& otherC) const{
      stack<monomial_t> stk;
    
      for(int i = 0; i < nMon; i++)
        compose(&mon[i], otherA, otherB, otherC, &stk);
    
      return polynomialFromStack(stk);
    }
    
    void Polynomial::operator=(const Polynomial& other){
      if(mon)
        delete[] mon;
    
      nMon = other.nMon;
      mon  = copyMonomial(other.mon, nMon);
    }
    
    string Polynomial::
    toString(const Polynomial::monomial_t* mon, const bool isAbs){
      stringstream stream;
      const bool minusOne    = mon->coef == -1.0;
      const bool notUnitCoef = mon->coef !=  1.0 && !minusOne;
    
      // If we have a constant term
      if(!mon->power[0] && !mon->power[1] && !mon->power[2]){
        stream << mon->coef;
        return stream.str();
      }
    
      // If we're here, we do not have a constant term
    
      // If we have a coefficient of '1', we don't display it
      if(notUnitCoef && isAbs)
        stream << abs(mon->coef);
    
      if(notUnitCoef && !isAbs)
        stream << mon->coef;
    
      if(minusOne && !isAbs)
        stream << "-";
    
      // We look for each power
      bool notOnce = false;
      for(int i = 0; i < 3; i++){
        // If we have a non zero power, we display it
        if(mon->power[i]){
          if(notUnitCoef || notOnce)
            stream << " * ";
    
          stream << coefName[i];
    
          if(mon->power[i] != 1)
            stream << "^" << mon->power[i];
    
          notOnce = true;
        }
      }
    
      return stream.str();
    }
    
    string Polynomial::toString(void) const{
      stringstream stream;
      bool isAbs = false;
    
      stream << toString(&mon[0], isAbs);
    
      for(int i = 1; i < nMon; i++){
        if(mon[i].coef < 0.0){
          stream << " - ";
          isAbs = true;
        }
    
        else{
          stream << " + ";
          isAbs = false;
        }
    
        stream << toString(&mon[i], isAbs);
      }
    
      return stream.str();
    }
    
    bool Polynomial::isSmaller(const Polynomial::monomial_t* a,
                               const Polynomial::monomial_t* b){
      // GRevLex order:
      // http://www.math.uiuc.edu/Macaulay2/doc/Macaulay2-1.4/share/doc/Macaulay2/Macaulay2Doc/html/___G__Rev__Lex.html
    
      int dif[3];
      int last = 0;
    
      if(isSmallerPower(a, b))
        return true;
    
      if(isEqualPower(a, b)){
        for(int i = 0, j = 2; i < 3; i++, j--)
          dif[i] = b->power[j] - a->power[j];
    
        for(int i = 0; i < 3; i++)
          if(dif[i])
            last = dif[i];
    
        if(last < 0)
          return true;
    
        else
          return false;
      }
    
      return false;
    }
    
    void Polynomial::sort(monomial_t* mon, int size){
      for(int i = 0; i < size; i++)
        for(int j = i; j < size; j++)
          if(isSmaller(&mon[j], &mon[i]))
            swap(mon, i, j);
    }
    
    void Polynomial::swap(monomial_t* mon, int i, int j){
      monomial_t tmp = mon[i];
      mon[i] = mon[j];
      mon[j] = tmp;
    }
    
    
    int Polynomial::mergeMon(monomial_t* sourceA, int sizeA,
                             monomial_t* sourceB, int sizeB,
                             monomial_t** dest){
      stack<monomial_t> s;
      monomial_t tmp;
    
      int i = 0;
      int j = 0;
      int N = 0;
    
      while(i < sizeA && j < sizeB){
        if(sourceA[i].coef == 0.0)
          i++;
    
        else if(sourceB[j].coef == 0.0)
          j++;
    
        else if(isEqual(&sourceA[i], &sourceB[j])){
          tmp       = sourceA[i];
          tmp.coef += sourceB[j].coef;
    
          if(tmp.coef != 0.0){
            s.push(tmp);
            N++;
          }
    
          i++;
          j++;
        }
    
        else if(isSmaller(&sourceA[i], &sourceB[j])){
          s.push(sourceA[i]);
          i++;
          N++;
        }
    
        else{
          s.push(sourceB[j]);
          j++;
          N++;
        }
      }
    
      while(i == sizeA && j < sizeB){
        s.push(sourceB[j]);
        j++;
        N++;
      }
    
      while(i < sizeA && j == sizeB){
        s.push(sourceA[i]);
        i++;
        N++;
      }
    
      if(!N){
        *dest = zeroPolynomial();
        N++;
      }
    
      else{
        *dest = new monomial_t[N];
    
        for(int k = N - 1; k >= 0; k--){
          (*dest)[k] = s.top();
          s.pop();
        }
      }
    
      return N;
    }
    
    int Polynomial::mult(const monomial_t* sourceA, int sizeA,
                         const monomial_t* sourceB, int sizeB,
                         monomial_t** dest){
    
      const monomial_t* a; // smaller polynomial
      const monomial_t* b; // bigger polynomial
      int nDist;
      int size;
    
      if(sizeA < sizeB){
        a     = sourceA;
        b     = sourceB;
        nDist = sizeA;
        size  = sizeB;
      }
    
      else{
        a     = sourceB;
        b     = sourceA;
        nDist = sizeB;
        size  = sizeA;
      }
    
      // Check if zero //
      if(a[0].coef == 0 || b[0].coef == 0){
        *dest = zeroPolynomial();
        return 1;
      }
    
      // Distrubute all monomials //
      monomial_t** dist = new monomial_t*[nDist];
    
      for(int i = 0; i < nDist; i++){
        dist[i] = copyMonomial(b, size);
    
        distribute(dist[i], size, &a[i]);
      }
    
      // Merge //
      int finalSize    = size;
      int nDistMinus   = nDist - 1;
      monomial_t** tmp = new monomial_t*[nDistMinus]; // Temp array for all dist[0];
    
      for(int i = 1, j = 0; i < nDist; i++, j++){
        tmp[j] = dist[0];
    
        finalSize = mergeMon(dist[0], finalSize,
                             dist[i], size,
                             &dist[0]);
      }
    
      // Keep distributed polynomial //
      *dest = dist[0];
    
      // Free Temporary Resources and Return //
      for(int i = 1, j = 0; i < nDist; i++, j++){
        delete[] dist[i];
        delete[] tmp[j];
      }
    
      delete[] dist;
      delete[] tmp;
    
      return finalSize;
    }
    
    void Polynomial::mult(monomial_t* source, int size, double alpha){
      for(int i = 0; i < size; i++)
        source[i].coef *= alpha;
    }
    
    
    void Polynomial::distribute(monomial_t* src, int size, const monomial_t* m){
      for(int i = 0; i < size; i++){
        src[i].coef *= m->coef;
    
        src[i].power[0] += m->power[0];
        src[i].power[1] += m->power[1];
        src[i].power[2] += m->power[2];
      }
    }
    
    void Polynomial::compose(const monomial_t* src,
                             Polynomial comp,
                             std::stack<monomial_t>* stk){
    
      comp.power(src->power[0]);
      comp.mul(src->coef);
    
      const int size = comp.nMon;
    
      for(int i = 0; i < size; i++){
        if(comp.mon[i].coef != 0){
    
          comp.mon[i].power[1] += src->power[1];
          comp.mon[i].power[2] += src->power[2];
    
          stk->push(comp.mon[i]);
        }
      }
    }
    
    void Polynomial::compose(const monomial_t* src,
                             Polynomial compA, Polynomial compB,
                             std::stack<monomial_t>* stk){
    
      compA.power(src->power[0]);
      compB.power(src->power[1]);
    
      compA.mul(compB);
      compA.mul(src->coef);
    
      const int size = compA.nMon;
    
      for(int i = 0; i < size; i++){
        if(compA.mon[i].coef != 0){
    
          compA.mon[i].power[2] += src->power[2];
    
          stk->push(compA.mon[i]);
        }
      }
    }
    
    void Polynomial::compose(const monomial_t* src,
                             Polynomial compA, Polynomial compB, Polynomial compC,
                             std::stack<monomial_t>* stk){
    
      compA.power(src->power[0]);
      compB.power(src->power[1]);
      compC.power(src->power[2]);
    
      compA.mul(compB);
      compA.mul(compC);
      compA.mul(src->coef);
    
      const int size = compA.nMon;
    
      for(int i = 0; i < size; i++){
        if(compA.mon[i].coef != 0)
          stk->push(compA.mon[i]);
      }
    }
    
    Polynomial Polynomial::
    polynomialFromStack(std::stack<Polynomial::monomial_t>& stk){
      Polynomial  newP;
      monomial_t* tmp;
      monomial_t* newMon;
      int         newNMon;
    
      if(!stk.size()){
        newMon  = zeroPolynomial();
        newNMon = 1;
      }
    
      else{
        newMon    = new monomial_t[1];
        newMon[0] = stk.top();
        newNMon   = 1;
        stk.pop();
    
        while(!stk.empty()){
          tmp     = newMon;
          newNMon = mergeMon(newMon, newNMon, &stk.top(), 1, &newMon);
          stk.pop();
    
          delete[] tmp;
        }
      }
    
      newP.nMon = newNMon;
      newP.mon  = newMon;
    
      return newP;
    }
    
    Polynomial::monomial_t* Polynomial::copyMonomial(const monomial_t* src,
                                                     int size){
      monomial_t* dest = new monomial_t[size];
    
      for(int i = 0; i < size; i++){
        dest[i].coef     = src[i].coef;
        dest[i].power[0] = src[i].power[0];
        dest[i].power[1] = src[i].power[1];
        dest[i].power[2] = src[i].power[2];
      }
    
      return dest;
    }
    
    Polynomial::monomial_t* Polynomial::zeroPolynomial(void){
      monomial_t* zero = new monomial_t[1];
    
      zero->coef     = 0;
      zero->power[0] = 0;
      zero->power[1] = 0;
      zero->power[2] = 0;
    
      return zero;
    }
    
    Polynomial::monomial_t* Polynomial::unitPolynomial(void){
      monomial_t* unit = new monomial_t[1];
    
      unit->coef     = 1;
      unit->power[0] = 0;
      unit->power[1] = 0;
      unit->power[2] = 0;
    
      return unit;
    }