Changeset 5797

Timestamp:

09/14/10 08:49:27 (15 years ago)

Author:

Mathieu Morlighem

Message:

Added Gauss points for quads

Location:

issm/trunk/src/c/objects/Gauss

Files:

• GaussPenta.cpp (modified) (3 diffs)
• GaussPenta.h (modified) (1 diff)

issm/trunk/src/c/objects/Gauss/GaussPenta.cpp ¶

@@ -92 +92 @@
         GaussLegendreLinear(&seg_coords,&seg_weights,numgauss);
 
-        /*Allocate GaussTria fields*/
+        /*Allocate GaussPenta fields*/
         coords1=(double*)xmalloc(numgauss*sizeof(double));
         coords2=(double*)xmalloc(numgauss*sizeof(double));
@@ -121 +121 @@
         }
         else{
-                ISSMERROR("Tria not supported yet");
+                ISSMERROR("Penta not supported yet");
         }
 
@@ -152 +152 @@
         else{
                 ISSMERROR("Tria not supported yet");
+        }
+
+}
+/*}}}*/
+/*FUNCTION GaussPenta::GaussPenta(int index1, int index2, int index3, int index4,int order_horiz,int order_vert){{{1*/
+GaussPenta::GaussPenta(int index1, int index2, int index3, int index4,int order_horiz,int order_vert){
+
+        /*Intermediaties*/
+        double *seg_horiz_coords  = NULL;
+        double *seg_horiz_weights = NULL;
+        double *seg_vert_coords   = NULL;
+        double *seg_vert_weights  = NULL;
+        int     i,j;
+
+        /*get the gauss points using the product of two line rules*/
+        GaussLegendreLinear(&seg_horiz_coords,&seg_horiz_weights,order_horiz);
+        GaussLegendreLinear(&seg_vert_coords, &seg_vert_weights, order_vert);
+
+        /*Allocate GaussPenta fields*/
+        numgauss=order_horiz*order_vert;
+        coords1=(double*)xmalloc(numgauss*sizeof(double));
+        coords2=(double*)xmalloc(numgauss*sizeof(double));
+        coords3=(double*)xmalloc(numgauss*sizeof(double));
+        coords4=(double*)xmalloc(numgauss*sizeof(double));
+        weights=(double*)xmalloc(numgauss*sizeof(double));
+
+        /*Quads: get the gauss points using the product of two line rules  */
+        if(index1==0 && index2==1 && index3==3 && index4==4){
+                for(i=0;i<order_horiz;i++){
+                        for(j=0;j<order_vert;j++){
+                                coords1[i*order_horiz+j]=  0.5*(1-seg_horiz_coords[i]);
+                                coords2[i*order_horiz+j]=1-0.5*(1-seg_horiz_coords[i]);
+                                coords3[i*order_horiz+j]=0.0;
+                                coords4[i*order_horiz+j]=seg_vert_coords[j];
+                                weights[i*order_horiz+j]=seg_horiz_weights[i]*seg_vert_weights[j];
+                        }
+                }
+        }
+        else if(index1==1 && index2==2 && index3==4 && index4==5){
+                for(i=0;i<order_horiz;i++){
+                        for(j=0;j<order_vert;j++){
+                                coords1[i*order_horiz+j]=0.0;
+                                coords2[i*order_horiz+j]=  0.5*(1-seg_horiz_coords[i]);
+                                coords3[i*order_horiz+j]=1-0.5*(1-seg_horiz_coords[i]);
+                                coords4[i*order_horiz+j]=seg_vert_coords[j];
+                                weights[i*order_horiz+j]=seg_horiz_weights[i]*seg_vert_weights[j];
+                        }
+                }
+        }
+        else if(index1==2 && index2==0 && index3==5 && index4==3){
+                for(i=0;i<order_horiz;i++){
+                        for(j=0;j<order_vert;j++){
+                                coords1[i*order_horiz+j]=1-0.5*(1-seg_horiz_coords[i]);
+                                coords2[i*order_horiz+j]=0.0;
+                                coords3[i*order_horiz+j]=  0.5*(1-seg_horiz_coords[i]);
+                                coords4[i*order_horiz+j]=seg_vert_coords[j];
+                                weights[i*order_horiz+j]=seg_horiz_weights[i]*seg_vert_weights[j];
+                        }
+                }
+        }
+        else{
+                ISSMERROR("Tria not supported yet (user provided indices %i %i %i %i)",index1,index2,index3,index4);
         }
 

issm/trunk/src/c/objects/Gauss/GaussPenta.h ¶

@@ -36 +36 @@
                 GaussPenta(int index1, int index2,int order);
                 GaussPenta(int index1, int index2, int index3, int order);
+                GaussPenta(int index1, int index2, int index3, int index4,int order_horiz,int order_vert);
                 ~GaussPenta();