Changeset 5738

Timestamp:

09/10/10 08:21:51 (15 years ago)

Author:

Mathieu Morlighem

Message:

Do not use beam in NumericalFlux

Location:

issm/trunk/src/c/objects

Files:

• Elements/Tria.cpp (modified) (1 diff)
• Elements/TriaRef.cpp (modified) (1 diff)
• Elements/TriaRef.h (modified) (1 diff)
• Loads/Numericalflux.cpp (modified) (8 diffs)

issm/trunk/src/c/objects/Elements/Tria.cpp

@@ -3884 +3884 @@
         }
 
-        /*Add Ke to global matrix Kgg: */
         MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES);
                 

issm/trunk/src/c/objects/Elements/TriaRef.cpp

@@ -118 +118 @@
                 *(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh3[0][i]; 
         }
+}
+/*}}}*/
+/*FUNCTION TriaRef::GetSegmentBFlux{{{1*/
+void TriaRef::GetSegmentBFlux(double* B,GaussTria* gauss, int index1,int index2){
+        /*Compute B  matrix. B=[phi1 phi2 -phi3 -phi4]
+         *
+         * and phi1=phi3 phi2=phi4
+         *
+         * We assume B has been allocated already, of size: 1x4
+         */
+
+        const int numgrids=3;
+        double l1l3[numgrids];
+
+        GetNodalFunctions(&l1l3[0],gauss);
+
+        B[0] = +l1l3[index1];
+        B[1] = +l1l3[index2];
+        B[2] = -l1l3[index1];
+        B[3] = -l1l3[index2];
+}
+/*}}}*/
+/*FUNCTION TriaRef::GetSegmentBprimeFlux{{{1*/
+void TriaRef::GetSegmentBprimeFlux(double* Bprime,GaussTria* gauss, int index1,int index2){
+        /*Compute Bprime  matrix. Bprime=[phi1 phi2 phi3 phi4]
+         *
+         * and phi1=phi3 phi2=phi4
+         *
+         * We assume Bprime has been allocated already, of size: 1x4
+         */
+
+        const int numgrids=3;
+        double l1l3[numgrids];
+
+        GetNodalFunctions(&l1l3[0],gauss);
+
+        Bprime[0] = l1l3[index1];
+        Bprime[1] = l1l3[index2];
+        Bprime[2] = l1l3[index1];
+        Bprime[3] = l1l3[index2];
 }
 /*}}}*/

issm/trunk/src/c/objects/Elements/TriaRef.h

@@ -48 +48 @@
                 void GetNodalFunctions(double* l1l2l3,GaussTria* gauss);
                 void GetSegmentNodalFunctions(double* l1l2l3,GaussTria* gauss, int index1,int index2);
+                void GetSegmentBFlux(double* B,GaussTria* gauss, int index1,int index2);
+                void GetSegmentBprimeFlux(double* Bprime,GaussTria* gauss, int index1,int index2);
                 void GetNodalFunctionsDerivatives(double* l1l2l3,double* xyz_list, double* gauss);
                 void GetNodalFunctionsDerivatives(double* l1l2l3,double* xyz_list, GaussTria* gauss);

issm/trunk/src/c/objects/Loads/Numericalflux.cpp

@@ -3 +3 @@
  */
 
+/*Headers:*/
+/*{{{1*/
 #ifdef HAVE_CONFIG_H
-        #include "config.h"
+#include "config.h"
 #else
 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
@@ -15 +17 @@
 #include "../../include/include.h"
 #include "../objects.h"
-
-extern int my_rank;
+/*}}}*/ 
+
+/*Load macros*/
+#define NUMVERTICES 4
+#define NDOF1 1
 
 /*Numericalflux constructors and destructor*/
@@ -394 +399 @@
 void  Numericalflux::CreateKMatrixInternal(Mat Kgg){
 
-        /* local declarations */
-        int             i,j,ig;
-        int analysis_type;
-
-        /* 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=NULL;
-
-        /* gaussian points: */
-        int     num_gauss;
-        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,vy;
-        double UdotN;
-        double dt;
-        int    found;
-        int    dim;
-
-        /*dynamic objects pointed to by hooks: */
-        Input  *vxaverage_input = NULL;
-        Input  *vyaverage_input = NULL;
-
-        /*Retrieve parameters: */
+        /* constants*/
+        const int numdof=NDOF1*NUMVERTICES;
+
+        /* Intermediaries*/
+        int        i,j,ig,index1,index2,analysis_type;
+        double     DL1,DL2,Jdet,dt,vx,vy,UdotN;
+        double     xyz_list[NUMVERTICES][3];
+        double     normal[2];
+        double     B[numdof];
+        double     Bprime[numdof];
+        double     Ke_gg1[numdof][numdof];
+        double     Ke_gg2[numdof][numdof];
+        double     Ke_gg[numdof][numdof] = {0.0};
+        int       *doflist = NULL;
+        GaussTria *gauss;
+
+        GetVerticesCoordinates(&xyz_list[0][0], nodes,NUMVERTICES);
+        GetDofList(&doflist);
+
+        /*Retrieve all inputs and parameters we will be needing: */
+        Tria*  tria=(Tria*)element;
+        Input* vxaverage_input=tria->inputs->GetInput(VxEnum);
+        Input* vyaverage_input=tria->inputs->GetInput(VyEnum);
         this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
-        parameters->FindParam(&dim,DimEnum);
-
-        /*recover objects from hooks: */
-        Tria* tria=(Tria*)element;
-
-        /*recover parameters: */
-        if (analysis_type==PrognosticAnalysisEnum){
-                parameters->FindParam(&dt,DtEnum);
-        }
-        else if (analysis_type==BalancedthicknessAnalysisEnum || analysis_type==AdjointBalancedthicknessAnalysisEnum){
-                /*No transient term is involved*/
-                dt=1;
-        }
-        else{
-                ISSMERROR("analysis_type %s not supported yet",EnumToString(analysis_type));
-        }
-
-        /* Get node coordinates, dof list and normal vector: */
-        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
-        GetDofList(&doflist);
         GetNormal(&normal[0],xyz_list);
-
-        /*Retrieve all inputs we will be needing: */
-        if(dim==2){
-                vxaverage_input=tria->inputs->GetInput(VxEnum);
-                vyaverage_input=tria->inputs->GetInput(VyEnum);
-        }
-
-        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
-        num_gauss=2;
-        gaussSegment(&gauss_coords, &gauss_weights,num_gauss);
+        switch(analysis_type){
+                case PrognosticAnalysisEnum:
+                        parameters->FindParam(&dt,DtEnum); break;
+                case BalancedthicknessAnalysisEnum: case AdjointBalancedthicknessAnalysisEnum:
+                        dt=1; break;/*No transient term is involved*/
+                default:
+                        ISSMERROR("analysis_type %s not supported yet",EnumToString(analysis_type));
+        }
 
         /* Start  looping on the number of gaussian points: */
-        for (ig=0; ig<num_gauss; ig++){
-
-                /*Pick up the gaussian point: */
-                gauss_weight=gauss_weights[ig];
-                gauss_coord =gauss_coords[ig];
-
-                /* Get Jacobian determinant: */
-                GetJacobianDeterminant(&Jdet,xyz_list,gauss_coord);
-
-                //Get vx, vy and v.n
-                this->GetParameterValue(&vx,gauss_coord,vxaverage_input);
-                this->GetParameterValue(&vy,gauss_coord,vyaverage_input);
-                
+        index1=tria->GetNodeIndex(nodes[0]);
+        index2=tria->GetNodeIndex(nodes[1]);
+        gauss=new GaussTria(index1,index2,2);
+        for(ig=gauss->begin();ig<gauss->end();ig++){
+
+                gauss->GaussPoint(ig);
+
+                tria->GetSegmentBFlux(&B[0],gauss,index1,index2);
+                tria->GetSegmentBprimeFlux(&Bprime[0],gauss,index1,index2);
+
+                vxaverage_input->GetParameterValue(&vx,gauss);
+                vyaverage_input->GetParameterValue(&vy,gauss);
                 UdotN=vx*normal[0]+vy*normal[1];
-        //      if (fabs(UdotN)<1.0e-9 && analysis_type==BalancedthicknessAnalysisEnum) printf("Edge number %i has a flux very small (u.n = %g ), which could lead to unaccurate results\n",id,UdotN);
-
-                /*Get L and B: */
-                GetL(&L[0],gauss_coord,numdof);
-                GetB(&B[0],gauss_coord);
-
-                /*Compute DLs*/
-                DL1=gauss_weight*Jdet*dt*UdotN/2;
-                DL2=gauss_weight*Jdet*dt*fabs(UdotN)/2;
-
-                /*  Do the triple products: */
+                tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
+                DL1=gauss->weight*Jdet*dt*UdotN/2;
+                DL2=gauss->weight*Jdet*dt*fabs(UdotN)/2;
+
                 TripleMultiply(&B[0],1,numdof,1,
                                         &DL1,1,1,0,
-                                        &L[0],1,numdof,0,
+                                        &Bprime[0],1,numdof,0,
                                         &Ke_gg1[0][0],0);
                 TripleMultiply(&B[0],1,numdof,1,
@@ -502 +460 @@
                                         &Ke_gg2[0][0],0);
 
-                /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
-                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg1[i][j];
-                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg2[i][j];
-
-        } // for (ig=0; ig<num_gauss; ig++)
-
-        /*Add Ke_gg to global matrix Kgg: */
+                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg1[i][j];
+                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg2[i][j];
+        }
+
         MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
-
-        xfree((void**)&gauss_coords);
-        xfree((void**)&gauss_weights);
         
         /*Free ressources:*/
+        delete gauss;
         xfree((void**)&doflist);
 
@@ -528 +481 @@
         /* node data: */
         const int numgrids=2;
-        const int NDOF1=1;
         const int numdof=NDOF1*numgrids;
         double    xyz_list[numgrids][3];
@@ -555 +507 @@
         double dt;
         int    dofs[1]={0};
-        int    found;
         int    dim;
 
@@ -668 +619 @@
         /* node data: */
         const int numgrids=2;
-        const int NDOF1=1;
         const int numdof=NDOF1*numgrids;
         double    xyz_list[numgrids][3];
@@ -695 +645 @@
         double dt;
         int    dofs[1]={0};
-        int    found;
         int    dim;