Index: /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp (revision 3370) +++ /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp (revision 3371) @@ -191,7 +191,7 @@ /*vertices offsets: */ - tria_g[0]=(int)*(iomodel->elements+elements_width*i+0); - tria_g[1]=(int)*(iomodel->elements+elements_width*i+1); - tria_g[2]=(int)*(iomodel->elements+elements_width*i+2); + tria_g[0]=3*i+1; + tria_g[1]=3*i+2; + tria_g[2]=3*i+3; /*thickness,surface and bed:*/ @@ -333,4 +333,5 @@ //Get id of the vertex on which the current node is located k=(int)*(iomodel->elements+elements_width*i+j)-1; //(Matlab to C indexing) + ISSMASSERT(k>=0 && knumberofnodes); //Get node properties Index: /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateLoadsPrognostic2.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateLoadsPrognostic2.cpp (revision 3370) +++ /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateLoadsPrognostic2.cpp (revision 3371) @@ -70,6 +70,6 @@ pos1=pos2=UNDEF; for(j=0;j<3;j++){ - if (iomodel->elements[3*(int)e1+j]==i1) pos1=j; - if (iomodel->elements[3*(int)e2+j]==i1) pos2=j; + if (iomodel->elements[3*(int)e1+j]==i1) pos1=j+1; + if (iomodel->elements[3*(int)e2+j]==i1) pos2=j+1; } ISSMASSERT(pos1!=UNDEF && pos2!=UNDEF); @@ -77,8 +77,8 @@ /*3: We have the id of the elements and the position of the vertices in the index * we can compute their dofs!*/ - numericalflux_node_ids[0]=3*(int)e1+pos1+1; //id is in matlab index - numericalflux_node_ids[1]=3*(int)e1+(pos1+1)%3+1; - numericalflux_node_ids[2]=3*(int)e2+pos2+1; - numericalflux_node_ids[3]=3*(int)e2+(pos2-1)%3+1; + numericalflux_node_ids[0]=3*(int)e1+pos1; //ex: 1 2 3 + numericalflux_node_ids[1]=3*(int)e1+(pos1%3)+1; //ex: 2 3 1 + numericalflux_node_ids[2]=3*(int)e2+pos2; //ex: 1 2 3 + numericalflux_node_ids[3]=3*(int)e2+((pos2+1)%3)+1; //ex: 3 1 2 } else{ @@ -94,5 +94,5 @@ pos1==UNDEF; for(j=0;j<3;j++){ - if (iomodel->elements[3*(int)e1+j]==i1) pos1=j; + if (iomodel->elements[3*(int)e1+j]==i1) pos1=j+1; } ISSMASSERT(pos1!=UNDEF); @@ -100,6 +100,6 @@ /*3: We have the id of the elements and the position of the vertices in the index * we can compute their dofs!*/ - numericalflux_node_ids[0]=3*(int)e1+pos1+1; //id is in matlab index - numericalflux_node_ids[1]=3*(int)e1+(pos1+1)%3+1; + numericalflux_node_ids[0]=3*(int)e1+pos1; + numericalflux_node_ids[1]=3*(int)e1+(pos1%3)+1; } @@ -107,5 +107,4 @@ loads->AddObject(numericalflux); - } Index: /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateParametersPrognostic2.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateParametersPrognostic2.cpp (revision 3370) +++ /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateParametersPrognostic2.cpp (revision 3371) @@ -74,5 +74,5 @@ IoModelFetchData(&thickness,NULL,NULL,iomodel_handle,"thickness"); h_g=(double*)xcalloc(iomodel->numberofelements*3,sizeof(double)); - if(thickness) for(i=0;i<3*iomodel->numberofelements;i++) h_g[i]=thickness[part[i]]/iomodel->yts; + if(thickness) for(i=0;i<3*iomodel->numberofelements;i++) h_g[i]=thickness[part[i]]; count++; @@ -85,5 +85,5 @@ IoModelFetchData(&accumulation,NULL,NULL,iomodel_handle,"accumulation"); a_g=(double*)xcalloc(iomodel->numberofelements*3,sizeof(double)); - if(accumulation) for(i=0;i<3*iomodel->numberofelements;i++) a_g[i]=accumulation[part[i]]/iomodel->yts; + if(accumulation) for(i=0;i<3*iomodel->numberofelements;i++) a_g[i]=accumulation[part[i]]; count++; @@ -96,5 +96,5 @@ IoModelFetchData(&melting,NULL,NULL,iomodel_handle,"melting"); m_g=(double*)xcalloc(iomodel->numberofelements*3,sizeof(double)); - if(melting) for(i=0;i<3*iomodel->numberofelements;i++) m_g[i]=melting[part[i]]/iomodel->yts; + if(melting) for(i=0;i<3*iomodel->numberofelements;i++) m_g[i]=melting[part[i]]; count++; Index: /issm/trunk/src/c/objects/Numericalflux.cpp =================================================================== --- /issm/trunk/src/c/objects/Numericalflux.cpp (revision 3370) +++ /issm/trunk/src/c/objects/Numericalflux.cpp (revision 3371) @@ -238,4 +238,20 @@ } /*}}}*/ +/*FUNCTION Numericalflux::GetB {{{1*/ +void Numericalflux::GetB(double* B, double gauss_coord){ + + const int numgrids=4; + double l1l4[numgrids]; + + /*Get nodal functions: */ + GetNodalFunctions(&l1l4[0],gauss_coord); + + /*Build B: */ + B[0] = +l1l4[0]; + B[1] = +l1l4[1]; + B[2] = -l1l4[0]; + B[3] = -l1l4[1]; +} +/*}}}*/ /*FUNCTION Numericalflux::GetDofList{{{1*/ @@ -254,5 +270,5 @@ } } - else if (strcmp(type,"boundary")){ + else if (strcmp(type,"boundary")==0){ for(i=0;i<2;i++){ nodes[i]->GetDofList(&doflist_per_node[0],&numberofdofspernode); @@ -262,5 +278,5 @@ } } - else ISSMERROR("type not supported yet"); + else ISSMERROR(exprintf("type %s not supported yet",type)); /*Assign output pointers:*/ @@ -288,4 +304,25 @@ } /*}}}*/ +/*FUNCTION Numericalflux::GetL {{{1*/ +void Numericalflux::GetL(double* L, double gauss_coord, int numdof){ + + const int numgrids=4; + double l1l4[numgrids]; + + /*Check numdof*/ + ISSMASSERT(numdof==2 || numdof==4); + + /*Get nodal functions: */ + GetNodalFunctions(&l1l4[0],gauss_coord); + + /*Luild L: */ + L[0] = l1l4[0]; + L[1] = l1l4[1]; + if (numdof==4){ + L[2] = l1l4[2]; + L[3] = l1l4[3]; + } +} +/*}}}*/ /*FUNCTION Numericalflux::GetName {{{1*/ char* Numericalflux::GetName(void){ @@ -300,7 +337,44 @@ l1l4[0]=-0.5*gauss_coord+0.5; l1l4[1]=+0.5*gauss_coord+0.5; - l1l4[0]=-0.5*gauss_coord+0.5; - l1l4[1]=+0.5*gauss_coord+0.5; - + l1l4[2]=-0.5*gauss_coord+0.5; + l1l4[3]=+0.5*gauss_coord+0.5; + +} +/*}}}*/ +/*FUNCTION Numericalflux::GetNormal {{{1*/ +void Numericalflux:: GetNormal(double* normal,double xyz_list[4][3]){ + + /*Build unit outward pointing vector*/ + const int numgrids=4; + double vector[2]; + double norm; + + vector[0]=xyz_list[1][0] - xyz_list[0][0]; + vector[1]=xyz_list[1][1] - xyz_list[0][1]; + + norm=sqrt(pow(vector[0],2.0)+pow(vector[1],2.0)); + + normal[0]= + vector[1]/norm; + normal[1]= - vector[0]/norm; +} +/*}}}*/ +/*FUNCTION Numericalflux::GetParameterValue {{{1*/ +void Numericalflux::GetParameterValue(double* pp, double* plist, double gauss_coord){ + + /*From node values of parameter p (plist[0],plist[1],plist[2]), return parameter value at gaussian + * point specifie by gauss_l1l2l3: */ + + /*nodal functions: */ + double l1l4[4]; + + /*output: */ + double p; + + GetNodalFunctions(&l1l4[0],gauss_coord); + + p=l1l4[0]*plist[0]+l1l4[1]*plist[1]; + + /*Assign output pointers:*/ + *pp=p; } /*}}}*/ @@ -326,14 +400,220 @@ /*}}}*/ /*FUNCTION Numericalflux::PenaltyCreateKMatrixInternal {{{1*/ -void Numericalflux::PenaltyCreateKMatrixInternal(Mat Kgg,void* inputs,double kmax,int analysis_type,int sub_analysis_type){ - - ISSMERROR("not supported yet"); - +void Numericalflux::PenaltyCreateKMatrixInternal(Mat Kgg,void* vinputs,double kmax,int analysis_type,int sub_analysis_type){ + + /* local declarations */ + int i,j; + + /* node data: */ + const int numgrids=4; + const int NDOF1=1; + const int numdof=NDOF1*numgrids; + double xyz_list[numgrids][3]; + double normal[2]; + int doflist[numdof]; + int numberofdofspernode; + + /* gaussian points: */ + int num_gauss,ig; + double* gauss_coords =NULL; + double gauss_coord; + double* gauss_weights=NULL; + double gauss_weight; + + /* matrices: */ + double B[numgrids]; + double L[numgrids]; + double DL1,DL2; + double Ke_gg1[numdof][numdof]; + double Ke_gg2[numdof][numdof]; + double Ke_gg[numdof][numdof]={0.0}; + double Jdet; + + /*input parameters for structural analysis (diagnostic): */ + double vx_list[numgrids]={0.0}; + double vy_list[numgrids]={0.0}; + double vx,vy; + double UdotN; + double dt; + int dofs[1]={0}; + int found; + + ParameterInputs* inputs=NULL; + + /*recover pointers: */ + inputs=(ParameterInputs*)vinputs; + + /*recover extra inputs from users, at current convergence iteration: */ + found=inputs->Recover("vx_average",&vx_list[0],1,dofs,numgrids,(void**)nodes); + if(!found)ISSMERROR(" could not find vx_average in inputs!"); + found=inputs->Recover("vy_average",&vy_list[0],1,dofs,numgrids,(void**)nodes); + if(!found)ISSMERROR(" could not find vy_average in inputs!"); + found=inputs->Recover("dt",&dt); + if(!found)ISSMERROR(" could not find dt in inputs!"); + + /* Get node coordinates, dof list and normal vector: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + GetNormal(&normal[0],xyz_list); + + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ + num_gauss=2; + GaussSegment(&gauss_coords, &gauss_weights,num_gauss); + + /* Start looping on the number of gaussian points: */ + for (ig=0; igRecover("vx_average",&vx_list[0],1,dofs,numgrids,(void**)nodes); + if(!found)ISSMERROR(" could not find vx_average in inputs!"); + found=inputs->Recover("vy_average",&vy_list[0],1,dofs,numgrids,(void**)nodes); + if(!found)ISSMERROR(" could not find vy_average in inputs!"); + found=inputs->Recover("dt",&dt); + if(!found)ISSMERROR(" could not find dt in inputs!"); + + /* Get node coordinates, dof list and normal vector: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + GetNormal(&normal[0],xyz_list); + + /*Check wether it is an inflow or outflow BC*/ + mean_vx=(vx_list[0]+vx_list[1])/2.0; + mean_vy=(vy_list[0]+vy_list[1])/2.0; + UdotN=mean_vx*normal[0]+mean_vy*normal[1]; + if (UdotN<=0){ + /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/ + return; + } + + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ + num_gauss=2; + GaussSegment(&gauss_coords, &gauss_weights,num_gauss); + + /* Start looping on the number of gaussian points: */ + for (ig=0; igRecover("vx_average",&vx_list[0],1,dofs,numgrids,(void**)nodes); + if(!found)ISSMERROR(" could not find vx_average in inputs!"); + found=inputs->Recover("vy_average",&vy_list[0],1,dofs,numgrids,(void**)nodes); + if(!found)ISSMERROR(" could not find vy_average in inputs!"); + found=inputs->Recover("dt",&dt); + if(!found)ISSMERROR(" could not find dt in inputs!"); + + /* Get node coordinates, dof list and normal vector: */ + GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + GetNormal(&normal[0],xyz_list); + + /*Check wether it is an inflow or outflow BC*/ + mean_vx=(vx_list[0]+vx_list[1])/2.0; + mean_vy=(vy_list[0]+vy_list[1])/2.0; + UdotN=mean_vx*normal[0]+mean_vy*normal[1]; + if (UdotN>0){ + /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/ + return; + } + + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ + num_gauss=2; + GaussSegment(&gauss_coords, &gauss_weights,num_gauss); + + /* Start looping on the number of gaussian points: */ + for (ig=0; ig