Select Git revision
robustPredicates.cpp
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;
}