Skip to content
Snippets Groups Projects
Select Git revision
  • 0a789359ee0e0db7e3d6b287d95fefc4c533f73d
  • master default protected
  • overlaps_tags_and_distributed_export
  • overlaps_tags_and_distributed_export_rebased
  • relaying
  • alphashapes
  • patches-4.14
  • steplayer
  • bl
  • pluginMeshQuality
  • fixBugsAmaury
  • hierarchical-basis
  • new_export_boris
  • oras_vs_osm
  • reassign_partitions
  • distributed_fwi
  • rename-classes
  • fix/fortran-api-example-t4
  • robust_partitions
  • reducing_files
  • fix_overlaps
  • gmsh_4_14_0
  • gmsh_4_13_1
  • gmsh_4_13_0
  • gmsh_4_12_2
  • gmsh_4_12_1
  • gmsh_4_12_0
  • gmsh_4_11_1
  • 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
41 results

robustPredicates.cpp

Blame
  • RandomField.cpp 13.69 KiB
    //
    // C++ random field class
    //
    // Description: random field for composite material parameters
    //
    // Author:  <L. Wu>, (C) 2016
    //
    // Copyright: See COPYING file that comes with this distribution
    //
    //
    #include <stdio.h>
    #include <stdlib.h>
    #include <math.h>
    #include <fstream>
    #include <string.h>
    #include "RandomField.h"
    #include "matrix_operations.h"
    
    #define PI 3.14159265
    //****************************************************************************
    //******************  class Yarn ***********************
    
    Yarn_class::Yarn_class(const double OrigX, const double OrigY, const double OrigZ, const double a, const double b, const double L, const double l, const char* WirGeo, const char* Warp_Weft): _OrigX(OrigX), _OrigY(OrigY), _OrigZ(OrigZ), _a(a), _b(b), _L(L), _l(l), _WirGeo(WirGeo), _Warp_Weft(Warp_Weft) {
    
         _c = sqrt(_a*_a -_b*_b);
    }
    
    Yarn_class::Yarn_class(const Yarn_class &source): _OrigX(source._OrigX), _OrigY(source._OrigY), _OrigZ(source._OrigZ), _a(source._a), _b(source._b), _c(source._c), _L(source._L),      _l(source._l),_WirGeo(source._WirGeo), _Warp_Weft(source._Warp_Weft) {};
    
    Yarn_class& Yarn_class::operator=(const Yarn_class &source){
    
       _OrigX = source._OrigX;
       _OrigY = source._OrigY;
       _OrigZ = source._OrigZ;
    
       _a = source._a;
       _b = source._b;
       _c = source._c;
       _L = source._L;
       _l = source._l;
       _WirGeo = source._WirGeo;
       _Warp_Weft = source._Warp_Weft;
    }
    
    Yarn_class::~Yarn_class(){};
    
    
    //**************** Gaussian Point in Yarn ?  **************************
    double Yarn_class::WireFunc(const double x, const double b, const double L) const{
        double tmp;
        if( _WirGeo == "Logistic"){   
            tmp = 1.0/(1.0+exp(-L*x));
            return 2.0*b*tmp-b;
        }
        else if(_WirGeo == "Cos"){ 
            tmp = 2.0*PI/L;
            return b*cos(tmp*x);
        }
    }
    
    double Yarn_class::WireFunc_Tang(const double x, const double b, const double L)const {
        double tmp;
        if( _WirGeo == "Logistic"){   
            tmp = 1.0/(1.0+exp(-L*x));
            return 2.0*b*tmp*(1.0-tmp)*L;
        }
        else if(_WirGeo == "Cos"){ 
            tmp = 2.0*PI/L;
            return -b*sin(tmp*x)*tmp;
        }
    }
    
    void Yarn_class::Inside_Yarn(const SVector3 &GaussP, double& k, double* V_Rand)const {
        double x, y, z;
        double x0, y0, z0;
        double tangent;
        double error = 1e-6;
        double gamma = 0.01;
        unsigned int max_iters = 10000;
        double L,b;
        x = GaussP(0) - _OrigX; 
        y = GaussP(1) - _OrigY; 
        z = GaussP(2); 
    
        if( _WirGeo == "Logistic"){
            L=_l;
            b = -_OrigZ;
            if(_Warp_Weft == "Weft") {
                if( x <= _L/2.0) x = x-_L/4.0;
                else{
                    x = x-3.0*_L/4.0;
                    b = -b;
                }
            } 
            else if(_Warp_Weft == "Warp"){
                if( y <= _L/2.0) y = y -_L/4.0;
                else{
                    y = y-3.0*_L/4.0;
                    b = -b;
                }
            }
        }
        else if(_WirGeo == "Cos"){
            b = _OrigZ;
            L=_L;
        }
    
        if(_Warp_Weft == "Weft"){
            x0 = x;           
            x0 = gradient_descent(x0, error, gamma, max_iters, x, z, b, L);
            z0 = WireFunc(x0, b, L);
            k = y*y/(_a*_a) + ((z-z0)*(z-z0)+(x-x0)*(x-x0))/(_b*_b);
            if(k<=1.0){
                double tangent = WireFunc_Tang(x0, b, L);
                V_Rand[0] = -90.0;
                V_Rand[1] = -atan(1.0/tangent) * 180.0 / PI;
                V_Rand[2] = 0.0;
            }
        }  
    
        else if(_Warp_Weft == "Warp"){
            y0 = y;           
            y0 = gradient_descent(y0, error, gamma, max_iters, y, z, b, L);
            z0 = WireFunc(y0, b, L);
            k = x*x/(_a*_a) + ((z-z0)*(z-z0)+(y-y0)*(y-y0))/(_b*_b);
            if(k<=1.0){
                tangent = WireFunc_Tang(y0, b, L);
                V_Rand[0] = 0.0;
                V_Rand[1] = -atan(1.0/tangent) * 180.0 / PI;
                V_Rand[2] = 0.0;
            }  
        }   
    }
    
    
    
    //******************************
    double Yarn_class::dfx(const double x, const double xG, const double yG, const double b, const double L) const {
        
        if(_WirGeo == "Cos"){
            double tmp = 2.0*PI/L;
            return 2.0*(x-xG) - 2.0*tmp*b*(b*cos(tmp*x)-yG)*sin(tmp*x);
         }
        else if(_WirGeo == "Logistic"){
            double tmp = 1.0/(1.0+exp(-L*x));
            return 2.0*(x-xG) - 2.0*(2.0*b*tmp-b-yG)*2.0*b*tmp*(1.0-tmp)*L;
        }
        else{
            printf("The center line function is not implemented! \n");
            return 0;
        }
    }
    
    
    double Yarn_class::gradient_descent(double x0, const double error, const double gamma, const unsigned int max_iters, const double xG, const double yG, const double b, 
                           const double L) const{
        double c_error = error + 1;
        unsigned int iters = 0;
        double p_error;
        while(error < c_error && iters < max_iters) {
            p_error = x0;
            x0 -= dfx(p_error, xG, yG, b, L) * gamma;
    	    c_error = fabs(p_error-x0);
            //printf("\nc_error %f\n", c_error);
            iters++;
        }
        return x0;
    }
    
    //**************************
    
    Yarn_class* init_Yarn(const double OrigX, const double OrigY, const double OrigZ, const double a, const double b, const double L, const double l, const char* WirGeo, const char* Warp_Weft){
    	Yarn_class* _Yarn;
    
    	_Yarn = new Yarn_class(OrigX, OrigY, OrigZ, a, b, L, l, WirGeo, Warp_Weft);
    			
            return _Yarn;
    }
    
    
    
    // general class for random field
    
    Random_Fclass::Random_Fclass(const SVector3& Ldomain, const double OrigX, const double OrigY, const double OrigZ, const double dx,
                                 const double dy, const double dz): _Ldomain(Ldomain), _OrigX(OrigX), _OrigY(OrigY), _OrigZ(OrigZ), 
                                 _dx(dx), _dy(dy), _dz(dz){
    
         _nx = Ldomain.x()/_dx+2;
         _ny = Ldomain.y()/_dy+2;
         _nz = Ldomain.z()/_dz+2;
    }
    
    Random_Fclass::Random_Fclass(const Random_Fclass &source): _Ldomain(source._Ldomain), _OrigX(source._OrigX), _OrigY(source._OrigY),
                                 _OrigZ(source._OrigZ), _dx(source._nx), _dy(source._dy), _dz(source._dz), _nx(source._nx), 
                                 _ny(source._ny), _nz(source._nz) {};
    
    Random_Fclass& Random_Fclass::operator=(const Random_Fclass &source){
       _Ldomain = source._Ldomain;
       _OrigX = source._OrigX;
       _OrigY = source._OrigY;
       _OrigZ = source._OrigZ;
    
       _dx = source._dx;
       _dy = source._dy;
       _dz = source._dz;
       _nx = source._nx;
       _ny = source._ny;
       _nz = source._nz;
    }
    
    Random_Fclass::~Random_Fclass(){};
    
    
    //****************  Random Euler Angler  **************************
    Random_Feul::Random_Feul(const double* mean, const SVector3& Ldomain, const double OrigX, const double OrigY, const double OrigZ, 
                  const double dx, const double dy, const double dz, const char* Geo):Random_Fclass(Ldomain, OrigX, OrigY, OrigZ, dx, dy, dz){
    
          mallocvector(&_Eul_mean,3);
          _Eul_mean[0] = mean[0];
          _Eul_mean[1] = mean[1];
          _Eul_mean[2] = mean[2];
          if (Geo!= NULL ) _Geo = Geo;
          if(_Geo == "Cos" or _Geo == "Logistic"){
              double a = dx;
              double b = dy;
              double L = _Ldomain.x();
              double l = dz;
              double tmpX, tmpY;
              char* weft = "Weft";
              char* warp = "Warp";
              tmpY = _OrigY - L/2.0;
              Weft1 = init_Yarn(_OrigX, tmpY, _OrigZ, a, b, L, l, _Geo,  weft);
              Weft2 = init_Yarn(_OrigX, _OrigY, _OrigZ, a, b, L, l, _Geo,  weft);
              tmpY = _OrigY + L/2.0;
              Weft3 = init_Yarn(_OrigX, tmpY, _OrigZ, a, b, L, l, _Geo,  weft);
    
              tmpY = _OrigY - L/2.0;
              Warp1 = init_Yarn(_OrigX, tmpY, _OrigZ, a, b, L, l, _Geo, warp);
              tmpX = _OrigX + L/2.0;
              Warp2 = init_Yarn(tmpX, tmpY, _OrigZ, a, b, L, l, _Geo, warp);
              tmpX = _OrigX + L;
              Warp3 = init_Yarn(tmpX, tmpY, _OrigZ, a, b, L, l, _Geo,  warp);  
          }
          malloctens4(&_RF_eul,_nx,_ny,_nz,3);
    
    }
    
    Random_Feul::Random_Feul(const Random_Feul &source):Random_Fclass(source){
    
          mallocvector(&_Eul_mean,3);
          _Eul_mean[0] = source._Eul_mean[0];
          _Eul_mean[1] = source._Eul_mean[1];
          _Eul_mean[2] = source._Eul_mean[2];
          _Geo = source._Geo;
    
          malloctens4(&_RF_eul,_nx,_ny,_nz,3);
          for(int i=0;i<_nx;i++){
              for(int j=0;j<_ny;j++){
                  for(int k=0;k<_nz;k++){
                      for(int p=0;p<3;p++){
                          _RF_eul[i][j][k][p] = source._RF_eul[i][j][k][p];
                      }
                  }
              }
          }
         Weft1 = NULL;
         Weft2 = NULL;
         Weft3 = NULL;
    
         Warp1 = NULL;
         Warp2 = NULL;
         Warp3 = NULL;
         if(source.Weft1 != NULL){
             Weft1 =source.Weft1->clone();
             Weft2 =source.Weft2->clone();    
             Weft3 =source.Weft3->clone();
    
             Warp1 =source.Warp1->clone();
             Warp2 =source.Warp2->clone();
             Warp3 =source.Warp3->clone();
         }
    }
    
    
    Random_Feul& Random_Feul::operator=(const Random_Fclass &source){
         Random_Fclass::operator=(source);
         const Random_Feul* src = dynamic_cast<const Random_Feul*>(&source);
         if(src!=NULL)
         {
            _Eul_mean[0] = src->_Eul_mean[0];
            _Eul_mean[1] = src->_Eul_mean[1];
            _Eul_mean[2] = src->_Eul_mean[2];
            _Geo = src->_Geo;
    
            for(int i=0;i<_nx;i++){
              for(int j=0;j<_ny;j++){
                  for(int k=0;k<_nz;k++){
                      for(int p=0;p<3;p++){
                          _RF_eul[i][j][k][p] = src->_RF_eul[i][j][k][p];
                      }
                  }
              }
            }
         }
    
         Weft1 = NULL;
         Weft2 = NULL;
         Weft3 = NULL;
    
         Warp1 = NULL;
         Warp2 = NULL;
         Warp3 = NULL;
         if(src->Weft1 != NULL){
             Weft1 =src->Weft1->clone();
             Weft2 =src->Weft2->clone();    
             Weft3 =src->Weft3->clone();
    
             Warp1 =src->Warp1->clone();
             Warp2 =src->Warp2->clone();
             Warp3 =src->Warp3->clone();
         }
    
    }
    
      
    Random_Feul::~Random_Feul(){
        free(_Eul_mean);
        freetens4(_RF_eul,_nx,_ny,_nz);
        if(Weft1 != NULL) {
             delete Weft1;
             delete Weft2;
             delete Weft3;
             delete Warp1;
             delete Warp2;
             delete Warp3;
        }
    }
    
    
    void Random_Feul::RandomGen(const SVector3 &GaussP, double* V_Rand){
        
         if(strcmp(_Geo,"Cylinder")==1)
         {
             if(_Eul_mean[0] == 0.0){
                 V_Rand[0] = 0.0;
                 V_Rand[1] = 0.0;
                 V_Rand[2] = 0.0;
              }
    
             else if(_Eul_mean[0] != 0.0){
                 double R = sqrt(pow((GaussP(0)-_OrigX), 2.0)+pow((GaussP(1)-_OrigY), 2.0));
                 V_Rand[0] = acos((GaussP(0)-_OrigX)/R)*180.0/PI;
                 if((GaussP(1)-_OrigY) < 0.0){
                     V_Rand[0] = -V_Rand[0];
                 }
                 V_Rand[1] = _Eul_mean[0];
                 V_Rand[2] = 0.0;
              }
         }
         else if(strcmp(_Geo,"Cos")==1 or strcmp(_Geo,"Logistic")==1){
             double x, y;
             double L = _Ldomain.x();
             const Yarn_class*  TWeft, *TWarp;
             
             x = GaussP(0)-_OrigX;
             y = GaussP(1)-_OrigY;
             if(x < L/4.0) TWarp = Warp1;
             else if(x > 3.0*L/4.0) TWarp = Warp3;
             else TWarp = Warp2;
      
             if(y < -L/4.0) TWeft = Weft1;
             else if(y > L/4.0) TWeft = Weft3;
             else TWeft = Weft2;
    
             double k_Warp, k_Weft;
             TWeft->Inside_Yarn( GaussP, k_Weft, V_Rand);
             TWarp->Inside_Yarn( GaussP, k_Warp, _Eul_mean);         
             if(k_Warp < k_Weft){
                 V_Rand[0] = _Eul_mean[0];
                 V_Rand[1] = _Eul_mean[1];
                 V_Rand[2] = _Eul_mean[2];
             
             } 
         }
    }
    
    //************** Random Material Propety *********************
    Random_FProp::Random_FProp(const int Randnum, const SVector3& Ldomain, const double OrigX, const double OrigY, const double OrigZ, 
             const double dx, const double dy, const double dz, const char* RandProp):Random_Fclass(Ldomain, OrigX, OrigY, OrigZ, 
             dx, dy, dz), _Randnum(Randnum) {
    
          malloctens4(&_RF_Prop,_nx,_ny,_nz,_Randnum);
    
        // function to fill _RF_Vfi with data file RandProp to generate ??????????????
    }
    
    Random_FProp::Random_FProp(const Random_FProp &source): Random_Fclass(source), _Randnum(source._Randnum) {
    
          malloctens4(&_RF_Prop,_nx,_ny,_nz,_Randnum);
          for(int i=0;i<_nx;i++){
              for(int j=0;j<_ny;j++){
                  for(int k=0;k<_nz;k++){
                      for(int p=0;p<_Randnum;p++){
                          _RF_Prop[i][j][k][p] = source._RF_Prop[i][j][k][p];
                      }
                  }
              }
          }
    }
    
    Random_FProp& Random_FProp::operator=(const Random_Fclass &source){
         Random_Fclass::operator=(source);
         const Random_FProp* src = dynamic_cast<const Random_FProp*>(&source);
         if(src!=NULL)
         {
            _Randnum = src->_Randnum;
    
            for(int i=0;i<_nx;i++){
               for(int j=0;j<_ny;j++){
                  for(int k=0;k<_nz;k++){
                      for(int p=0;p<_Randnum;p++){
                          _RF_Prop[i][j][k][p] = src->_RF_Prop[i][j][k][p];
                      }
                  }
               }
            }
         }
    }
    
    Random_FProp::~Random_FProp(){
          freetens4(_RF_Prop,_nx,_ny,_nz);
    }
    
    void Random_FProp::RandomGen(const SVector3 &GaussP, double* V_Rand){
    
        double V = (rand()/(double)RAND_MAX);
        double _Vfi_mean = 0.5;
    
        V_Rand[0] = _Vfi_mean + (V-0.5)*0.05;
    }
    
    
    
    
    //**************************
    
    Random_Fclass* init_RandF(const double* mean, const SVector3& Ldomain, const double OrigX, const double OrigY, const double OrigZ,
                              const double dx, const double dy, const double dz, const char* Geo){
    	Random_Fclass* _RFclass;
    
            _RFclass = new Random_Feul(mean, Ldomain, OrigX, OrigY, OrigZ, dx, dy, dz, Geo);
    
            return _RFclass;
    }
    
    
    
    Random_Fclass* init_RandF(const int Randnum, const SVector3& Ldomain, const double OrigX, const double OrigY, const double OrigZ,
                              const double dx, const double dy, const double dz, const char*RandProp){
    	Random_Fclass* _RFclass;
    
    	_RFclass = new Random_FProp(Randnum, Ldomain, OrigX, OrigY, OrigZ, dx, dy, dz, RandProp);
    			
            return _RFclass;
    }