Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 4549) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 4550) @@ -32,4 +32,27 @@ } /*}}}*/ +/*FUNCTION Tria::Tria(int id, int index, IoModel* iomodel,int nummodels){{{1*/ +Tria::Tria(int tria_id, int index, IoModel* iomodel,int nummodels) + + :TriaRef(P1Enum) //P1Enum: interpolation type + ,TriaHook(nummodels,index+1,iomodel->numberofelements+1) //index+1: matice id, iomodel->numberofelements+1: matpar id + { //i is the element index + /*id: */ + this->id=tria_id; + + //this->parameters: we still can't point to it, it may not even exist. Configure will handle this. + this->parameters=NULL; + + /*intialize inputs and results: */ + this->inputs=new Inputs(); + this->results=new Results(); + + /*initialize pointers:*/ + this->nodes=NULL; + this->matice=NULL; + this->matpar=NULL; + +} +/*}}}*/ /*FUNCTION Tria::~Tria(){{{1*/ Tria::~Tria(){ @@ -37,27 +60,4 @@ delete results; this->parameters=NULL; -} -/*}}}*/ -/*FUNCTION Tria::Tria(int id, int index, IoModel* iomodel,int nummodels){{{1*/ -Tria::Tria(int tria_id, int index, IoModel* iomodel,int nummodels) - - :TriaRef(P1Enum) //P1Enum: interpolation type - ,TriaHook(nummodels,index+1,iomodel->numberofelements+1) //index+1: matice id, iomodel->numberofelements+1: matpar id - { //i is the element index - /*id: */ - this->id=tria_id; - - //this->parameters: we still can't point to it, it may not even exist. Configure will handle this. - this->parameters=NULL; - - /*intialize inputs and results: */ - this->inputs=new Inputs(); - this->results=new Results(); - - /*initialize pointers:*/ - this->nodes=NULL; - this->matice=NULL; - this->matpar=NULL; - } /*}}}*/ Index: /issm/trunk/src/c/objects/Materials/Matice.cpp =================================================================== --- /issm/trunk/src/c/objects/Materials/Matice.cpp (revision 4549) +++ /issm/trunk/src/c/objects/Materials/Matice.cpp (revision 4550) @@ -19,61 +19,69 @@ /*FUNCTION Matice::Matice(){{{1*/ Matice::Matice(){ + this->inputs=NULL; return; } /*}}}*/ -/*FUNCTION Matice::Matice(int in_mid,double in_B,double in_n){{{1*/ -Matice::Matice(int in_mid,double in_B,double in_n){ - this->Init(in_mid,in_B,in_n); -} -/*}}}*/ -/*FUNCTION Matice::Matice(int id, int i, IoModel* iomodel, int num_vertices){{{1*/ -Matice::Matice(int matice_mid,int i, IoModel* iomodel){ - - int j; - - /*needed for Init routine:*/ - double matice_B; - double matice_n; - int num_vertices; - - /*intermediary: */ - double B_avg; - - /*2d or 3d? */ +/*FUNCTION Matice::Matice(int id, int index, IoModel* iomodel, int num_vertices){{{1*/ +Matice::Matice(int matice_mid,int index, IoModel* iomodel){ + + /*Intermediaries:*/ + int i; + + /*Initialize id*/ + this->mid=matice_mid; + + /*Initialize inputs*/ + this->inputs=new Inputs(); + + /*if 2d*/ if(iomodel->dim==2){ - num_vertices=3; //tria elements - } + + /*Intermediaries*/ + const int num_vertices = 3; //Tria has 3 vertices + double nodeinputs[num_vertices]; + + /*Get B*/ + if (iomodel->rheology_B) { + for(i=0;irheology_B[int(iomodel->elements[num_vertices*index+i]-1)]; + this->inputs->AddInput(new TriaVertexInput(RheologyBEnum,nodeinputs)); + } + + /*Get n*/ + if (iomodel->rheology_n) { + for(i=0;irheology_n[int(iomodel->elements[num_vertices*index+i]-1)]; + this->inputs->AddInput(new TriaVertexInput(RheologyNEnum,nodeinputs)); + } + } + + /*if 3d*/ else if(iomodel->dim==3){ - num_vertices=6; //penta elements - } + + /*Intermediaries*/ + const int num_vertices = 6; //Penta has 6 vertices + double nodeinputs[num_vertices]; + + /*Get B*/ + if (iomodel->rheology_B) { + for(i=0;irheology_B[int(iomodel->elements[num_vertices*index+i]-1)]; + this->inputs->AddInput(new PentaVertexInput(RheologyBEnum,nodeinputs)); + } + + /*Get n*/ + if (iomodel->rheology_n) { + for(i=0;irheology_n[int(iomodel->elements[num_vertices*index+i]-1)]; + this->inputs->AddInput(new PentaVertexInput(RheologyNEnum,nodeinputs)); + } + } + + /*Else*/ else ISSMERROR(" Mesh type not supported yet!"); - /*Compute B on the element if provided*/ - if (iomodel->rheology_B){ - B_avg=0; - for(j=0;jrheology_B+((int)*(iomodel->elements+num_vertices*i+j)-1)); - } - B_avg=B_avg/num_vertices; - matice_B=B_avg; - } - else matice_B=UNDEF; - - if (iomodel->rheology_n) matice_n=(double)*(iomodel->rheology_n+i); - else matice_n=UNDEF; - - this->Init(matice_mid,matice_B,matice_n); } /*}}}*/ /*FUNCTION Matice::~Matice(){{{1*/ Matice::~Matice(){ + delete inputs; return; -} -/*}}}*/ -/*FUNCTION Matice::Init {{{1*/ -void Matice::Init(int in_mid,double in_B,double in_n){ - this->mid=in_mid; - this->B=in_B; - this->n=in_n; } /*}}}*/ @@ -85,7 +93,6 @@ printf("Matice:\n"); printf(" mid: %i\n",mid); - printf(" B: %g\n",B); - printf(" n: %g\n",n); - return; + printf(" inputs:\n"); + inputs->Echo(); } /*}}}*/ @@ -95,7 +102,6 @@ printf("Matice:\n"); printf(" mid: %i\n",mid); - printf(" B: %g\n",B); - printf(" n: %g\n",n); - return; + printf(" inputs:\n"); + inputs->DeepEcho(); } /*}}}*/ @@ -112,6 +118,9 @@ void Matice::Marshall(char** pmarshalled_dataset){ + /*Intermediaries*/ char* marshalled_dataset=NULL; int enum_type=0; + char* marshalled_inputs=NULL; + int marshalled_inputs_size; /*recover marshalled_dataset: */ @@ -126,9 +135,15 @@ /*marshall Matice data: */ memcpy(marshalled_dataset,&mid,sizeof(mid));marshalled_dataset+=sizeof(mid); - memcpy(marshalled_dataset,&B,sizeof(B));marshalled_dataset+=sizeof(B); - memcpy(marshalled_dataset,&n,sizeof(n));marshalled_dataset+=sizeof(n); + + /*Marshall inputs: */ + marshalled_inputs_size=inputs->MarshallSize(); + marshalled_inputs=inputs->Marshall(); + memcpy(marshalled_dataset,marshalled_inputs,marshalled_inputs_size*sizeof(char)); + marshalled_dataset+=marshalled_inputs_size; *pmarshalled_dataset=marshalled_dataset; - return; + + /*clean up and return*/ + xfree((void**)&marshalled_inputs); } /*}}}*/ @@ -136,5 +151,7 @@ int Matice::MarshallSize(){ - return sizeof(mid)+sizeof(B)+sizeof(n)+sizeof(int); //sizeof(int) for enum type + return sizeof(mid) + +inputs->MarshallSize() + +sizeof(int); //sizeof(int) for enum type } /*}}}*/ @@ -151,6 +168,7 @@ memcpy(&mid,marshalled_dataset,sizeof(mid));marshalled_dataset+=sizeof(mid); - memcpy(&B,marshalled_dataset,sizeof(B));marshalled_dataset+=sizeof(B); - memcpy(&n,marshalled_dataset,sizeof(n));marshalled_dataset+=sizeof(n); + + /*demarshall inputs: */ + inputs=(Inputs*)DataSetDemarshallRaw(&marshalled_dataset); /*return: */ @@ -168,5 +186,16 @@ /*FUNCTION Matice::copy {{{1*/ Object* Matice::copy() { - return new Matice(*this); + + /*Output*/ + Matice* matice=NULL; + + matice->mid=this->mid; + if(this->inputs){ + matice->inputs=(Inputs*)this->inputs->Copy(); + } + else{ + matice->inputs=new Inputs(); + } + return matice; } /*}}}*/ @@ -175,4 +204,9 @@ /*FUNCTION Matice::GetB {{{1*/ double Matice::GetB(){ + + /*Output*/ + double B; + + inputs->GetParameterAverage(&B,RheologyBEnum); return B; } @@ -180,4 +214,9 @@ /*FUNCTION Matice::GetN {{{1*/ double Matice::GetN(){ + + /*Output*/ + double n; + + inputs->GetParameterAverage(&n,RheologyNEnum); return n; } @@ -203,6 +242,11 @@ double exx,eyy,exy; - /*Intermediary value A and exponent e: */ + /*Intermediary: */ double A,e; + double B,n; + + /*Get B and n*/ + B=GetB(); + n=GetN(); if (n==1){ @@ -264,6 +308,11 @@ double exx,eyy,exy,exz,eyz; - /*Intermediary value A and exponent e: */ + /*Intermediaries: */ double A,e; + double B,n; + + /*Get B and n*/ + B=GetB(); + n=GetN(); if (n==1){ @@ -310,10 +359,9 @@ /*FUNCTION Matice::GetViscosity3dStokes {{{1*/ void Matice::GetViscosity3dStokes(double* pviscosity3d, double* epsilon){ - /*Return viscosity accounting for steady state power law creep [Thomas and MacAyeal, 1982]: * * 2*B * viscosity3d= ------------------------------------------------------------------- - * 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n] + * 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n] * * where mu is the viscotiy, B the flow law parameter , (u,v) the velocity @@ -330,9 +378,13 @@ double exx,eyy,exy,exz,eyz,ezz; - /*Intermediary value A and exponent e: */ + /*Intermediaries: */ double A,e; + double B,n; double eps0; + /*Get B and n*/ eps0=pow((double)10,(double)-27); + B=GetB(); + n=GetN(); if (n==1){ @@ -379,5 +431,4 @@ /*FUNCTION Matice::GetViscosityComplement {{{1*/ void Matice::GetViscosityComplement(double* pviscosity_complement, double* epsilon){ - /*Return viscosity accounting for steady state power law creep [Thomas and MacAyeal, 1982]: * @@ -401,4 +452,9 @@ /*Intermediary value A and exponent e: */ double A,e; + double B,n; + + /*Get B and n*/ + B=GetB(); + n=GetN(); if(epsilon){ @@ -432,9 +488,4 @@ } /*}}}*/ -/*FUNCTION Matice::SetB {{{1*/ -void Matice::SetB(double B_param){ - B=B_param; -} -/*}}}*/ /*FUNCTION Matice::InputUpdateFromVector(double* vector, int name, int type) {{{1*/ void Matice::InputUpdateFromVector(double* vector, int name, int type){ Index: /issm/trunk/src/c/objects/Materials/Matice.h =================================================================== --- /issm/trunk/src/c/objects/Materials/Matice.h (revision 4549) +++ /issm/trunk/src/c/objects/Materials/Matice.h (revision 4550) @@ -16,6 +16,5 @@ private: int mid; - double B; - double n; + Inputs* inputs; public: @@ -23,7 +22,5 @@ /*Matice constructors, destructors: {{{1*/ Matice(); - Matice(int mid,double B,double n); Matice(int mid,int i, IoModel* iomodel); - void Init(int mid,double B,double n); ~Matice(); /*}}}*/ @@ -49,5 +46,4 @@ /*}}}*/ /*Matice Numerics: {{{1*/ - void SetB(double B_param); void GetViscosity2d(double* pviscosity, double* pepsilon); void GetViscosity3d(double* pviscosity3d, double* pepsilon);