Changeset 25451


Ignore:
Timestamp:
08/24/20 09:22:09 (5 years ago)
Author:
tsantos
Message:

CHG: methods to compute discontinuous slope for surface elevation

Location:
issm/trunk-jpl/src/c
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp

    r25438 r25451  
    2020//#define FSANALYTICAL 10
    2121//#define LATERALFRICTION 1
     22//#define DISCSLOPE 1 //testing for SSA
    2223
    2324/*Model processing*/
    … …  
    17081709
    17091710        /* Start  looping on the number of gaussian points: */
     1711        #ifndef DISCSLOPE
    17101712        Gauss* gauss=element->NewGauss(2);
     1713        #else
     1714        Gauss* gauss=NULL;
     1715        Gauss* gauss_subelem=NULL;
     1716        bool partly_floating=false;
     1717        bool mainlyfloating=false;
     1718        int point1;
     1719        IssmDouble fraction1,fraction2;
     1720        IssmDouble phi=element->GetGroundedPortion(xyz_list);
     1721        if(phi>0.00000001 && phi<0.99999999) partly_floating=true;
     1722
     1723        int ig=-1;// needed for driving stress parameterization
     1724        if(partly_floating){
     1725                element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
     1726           gauss=element->NewGauss(point1,fraction1,fraction2,3); //considering the entire element
     1727                gauss_subelem=element->NewGauss(fraction1,fraction2,3);//gauss on each subelement
     1728        }
     1729        else{
     1730                gauss=element->NewGauss(2);//original
     1731        }
     1732        #endif
     1733
    17111734        while(gauss->next()){
    17121735
    … …  
    17161739                thickness_input->GetInputValue(&thickness,gauss); _assert_(thickness>0);
    17171740                surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
     1741                #ifdef DISCSLOPE
     1742                if(gauss_subelem && partly_floating){
     1743                        ig++;
     1744                        gauss_subelem->next();
     1745                        _assert_(std::abs(gauss_subelem->weight-gauss->weight)<0.0000001);
     1746                        /*Compute the discontinuous surface slope for this gauss point/subelement*/
     1747                        this->ComputeSurfaceSlope(&slope[0],gauss_subelem,gauss,point1,fraction1,fraction2,ig,dim,element);
     1748                }               
     1749                #endif
    17181750
    17191751                for(int i=0;i<numnodes;i++){
    … …  
    17301762        xDelete<IssmDouble>(basis);
    17311763        delete gauss;
     1764        #ifdef DISCSLOPE
     1765        if(gauss_subelem) delete gauss_subelem;
     1766        #endif
    17321767        return pe;
    17331768}/*}}}*/
    … …  
    20332068        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    20342069}/*}}}*/
     2070#ifdef DISCSLOPE
     2071void StressbalanceAnalysis::ComputeSurfaceSlope(IssmDouble* slope,Gauss* gauss_DG,Gauss* gauss_CG,int point1,IssmDouble fraction1,IssmDouble fraction2,int ig,int dim,Element* element){/*{{{*/
     2072
     2073        /*Compute the surface slope for each subelement, for a given integration point (gauss)*/
     2074        int numnodes=element->GetNumberOfNodes();
     2075        IssmDouble rho_ice=element->FindParam(MaterialsRhoIceEnum);
     2076        IssmDouble rho_water=element->FindParam(MaterialsRhoSeawaterEnum); //ocean
     2077        IssmDouble* H=xNew<IssmDouble>(numnodes);
     2078        IssmDouble* S=xNew<IssmDouble>(numnodes);
     2079        IssmDouble* S_subelem=xNew<IssmDouble>(numnodes);
     2080        IssmDouble H_f1,H_f2,S_f1,S_f2;
     2081
     2082        /*Get nodal deriviatives of the subelements related to the Gauss point*/
     2083        IssmDouble* dbasis_subelem=xNew<IssmDouble>(dim*numnodes);//CG basis for each subelement
     2084        this->NodalFunctionsDerivativesRGB(dbasis_subelem,gauss_DG,gauss_CG,point1,fraction1,fraction2,ig,dim,element);
     2085
     2086        /*Define thickness at the grounding line (on element edges)*/
     2087        element->GetInputListOnNodes(H,ThicknessEnum); // thickness on vertices
     2088        switch(point1){//{{{
     2089                case 0:
     2090                        H_f1=H[0]+(H[1]-H[0])*fraction1;
     2091                        H_f2=H[0]+(H[2]-H[0])*fraction2;
     2092                        break;
     2093                case 1:
     2094                        H_f1=H[1]+(H[2]-H[1])*fraction1;
     2095                        H_f2=H[1]+(H[0]-H[1])*fraction2;
     2096                        break;
     2097                case 2:
     2098                        H_f1=H[2]+(H[0]-H[2])*fraction1;
     2099                        H_f2=H[2]+(H[1]-H[2])*fraction2;
     2100                        break;
     2101                default:
     2102                        _error_("index "<<point1<<" not supported yet");
     2103        }//}}}
     2104
     2105        /*Define surface at the grounding (on element edges)*/
     2106        S_f1=H_f1*(1-rho_ice/rho_water);
     2107        S_f2=H_f2*(1-rho_ice/rho_water);
     2108        element->GetInputListOnNodes(S,SurfaceEnum); // surface on vertices
     2109
     2110        /*Define surface on the subelement vertices*/
     2111   if(ig<4){ // BLUE element itapopo only if order is = 3
     2112                switch(point1){ //{{{
     2113                        case 0:
     2114                                S_subelem[0]=S[0];
     2115                                S_subelem[1]=S_f1;
     2116                                S_subelem[2]=S_f2;
     2117                                break;
     2118                        case 1:
     2119                                S_subelem[0]=S[1];
     2120                                S_subelem[1]=S_f1;
     2121                                S_subelem[2]=S_f2;
     2122                                break;
     2123                        case 2:
     2124                                S_subelem[0]=S[2];
     2125                                S_subelem[1]=S_f1;
     2126                                S_subelem[2]=S_f2;
     2127                                break;
     2128                        default:
     2129                                _error_("index "<<point1<<" not supported yet");
     2130                }//}}} 
     2131        }
     2132        if(ig>3 && ig<8){ // GREEN element
     2133                switch(point1){ //{{{
     2134                        case 0:
     2135                                S_subelem[0]=S_f1;
     2136                                S_subelem[1]=S[2];
     2137                                S_subelem[2]=S_f2;
     2138                                break;
     2139                        case 1:
     2140                                S_subelem[0]=S_f1;
     2141                                S_subelem[1]=S[0];
     2142                                S_subelem[2]=S_f2;
     2143                                break;
     2144                        case 2:
     2145                                S_subelem[0]=S_f1;
     2146                                S_subelem[1]=S[1];
     2147                                S_subelem[2]=S_f2;
     2148                                break;
     2149                        default:
     2150                                _error_("index "<<point1<<" not supported yet");
     2151                }//}}} 
     2152        }
     2153        if(ig>7){ // RED element
     2154                switch(point1){ //{{{
     2155                        case 0:
     2156                                S_subelem[0]=S_f1;
     2157                                S_subelem[1]=S[1];
     2158                                S_subelem[2]=S[2];
     2159                                break;
     2160                        case 1:
     2161                                S_subelem[0]=S_f1;
     2162                                S_subelem[1]=S[2];
     2163                                S_subelem[2]=S[0];
     2164                                break;
     2165                        case 2:
     2166                                S_subelem[0]=S_f1;
     2167                                S_subelem[1]=S[0];
     2168                                S_subelem[2]=S[1];
     2169                                break;
     2170                        default:
     2171                                _error_("index "<<point1<<" not supported yet");
     2172                }//}}} 
     2173        }
     2174
     2175        //_printf_("\t"<<"H[1]\t"<<"H[2]\t"<<"H[3]\t"<<"S[1]\t"<<"S[2]\t"<<"S[3]\n");
     2176        //_printf_("\t"<<H[0]<<"\t"<<H[1]<<"\t"<<H[2]<<"\t"<<S[0]<<"\t"<<S[1]<<"\t"<<S[2]<<"\n");
     2177        //_printf_("\t H_f1 \t H_f2 \t S_f1 \t S_f2 \n");
     2178        //_printf_("\t"<< H_f1 <<"\t"<< H_f2 << "\t" << S_f1 << "\t" << S_f2<< "\n");
     2179        //_printf_("\t"<<"S_subelem[1]\t"<<"S_subelem[2]\t"<<"S_subelem[3]\n");
     2180        //_printf_("\t"<<S_subelem[0]<<"\t"<<S_subelem[1]<<"\t"<<S_subelem[2]<<"\n");
     2181
     2182        /*Compute slope*/
     2183        slope[0]=0;
     2184        slope[1]=0;
     2185        for(int i=0;i<numnodes;i++){
     2186                slope[0]+=S_subelem[i]*dbasis_subelem[numnodes*0+i]; //dSdx
     2187                slope[1]+=S_subelem[i]*dbasis_subelem[numnodes*1+i]; //dSdy
     2188        }
     2189
     2190        /*Clean up*/
     2191        xDelete<IssmDouble>(dbasis_subelem);
     2192        xDelete<IssmDouble>(H);
     2193        xDelete<IssmDouble>(S);
     2194        xDelete<IssmDouble>(S_subelem);
     2195
     2196}/*}}}*/
     2197void StressbalanceAnalysis::NodalFunctionsDerivativesRGB(IssmDouble* dbasis_subelem,Gauss* gauss_DG,Gauss* gauss_CG,int point1,IssmDouble fraction1,IssmDouble fraction2,int ig,int dim,Element* element){/*{{{*/
     2198       
     2199        /*Fetch number of nodes for this finite element*/
     2200   int numnodes = element->GetNumberOfNodes();
     2201   IssmDouble dbasis_ref[dim*numnodes];
     2202   IssmDouble Jinv[2][2];
     2203   //IssmDouble* dbasis_subelem = xNew<IssmDouble>(dim*numnodes);//CG basis for each subelement
     2204   IssmDouble* xyz_list_RGB = xNew<IssmDouble>(3*numnodes); // x,y,z per node
     2205   IssmDouble* xyz_list = NULL;
     2206
     2207   element->GetVerticesCoordinates(&xyz_list);
     2208   Tria* tria = xDynamicCast<Tria*>(element);
     2209
     2210   /*Get nodal functions derivatives*/
     2211   tria->GetNodalFunctionsDerivativesReference(dbasis_ref,gauss_DG,P1Enum); //gauss is not used here
     2212
     2213   if(ig<4){ // BLUE element itapopo only if order is = 3
     2214      switch(point1){//{{{
     2215         case 0:
     2216            /*Vertex 0 - local in the subelement*/
     2217            xyz_list_RGB[3*0+0]=xyz_list[3*0+0]; //x
     2218            xyz_list_RGB[3*0+1]=xyz_list[3*0+1]; //y
     2219            /*Vertex 1*/
     2220            xyz_list_RGB[3*1+0]=xyz_list[3*0+0]+(xyz_list[3*1+0]-xyz_list[3*0+0])*fraction1; //x
     2221            xyz_list_RGB[3*1+1]=xyz_list[3*0+1]+(xyz_list[3*1+1]-xyz_list[3*0+1])*fraction1; //y
     2222            /*Vertex 2*/
     2223            xyz_list_RGB[3*2+0]=xyz_list[3*0+0]+(xyz_list[3*2+0]-xyz_list[3*0+0])*fraction2; //x;
     2224            xyz_list_RGB[3*2+1]=xyz_list[3*0+1]+(xyz_list[3*2+1]-xyz_list[3*0+1])*fraction2; //y;
     2225            break;
     2226                         case 1:
     2227            /*Vertex 0 - local in the subelement*/
     2228            xyz_list_RGB[3*0+0]=xyz_list[3*1+0]; //x
     2229            xyz_list_RGB[3*0+1]=xyz_list[3*1+1]; //y
     2230            /*Vertex 1*/
     2231            xyz_list_RGB[3*1+0]=xyz_list[3*1+0]+(xyz_list[3*2+0]-xyz_list[3*1+0])*fraction1; //x
     2232            xyz_list_RGB[3*1+1]=xyz_list[3*1+1]+(xyz_list[3*2+1]-xyz_list[3*1+1])*fraction1; //y
     2233            /*Vertex 2*/
     2234            xyz_list_RGB[3*2+0]=xyz_list[3*1+0]+(xyz_list[3*0+0]-xyz_list[3*1+0])*fraction2; //x;
     2235            xyz_list_RGB[3*2+1]=xyz_list[3*1+1]+(xyz_list[3*0+1]-xyz_list[3*1+1])*fraction2; //y;
     2236            break;
     2237                        case 2:
     2238            /*Vertex 0 - local in the subelement*/
     2239            xyz_list_RGB[3*0+0]=xyz_list[3*2+0]; //x
     2240            xyz_list_RGB[3*0+1]=xyz_list[3*2+1]; //y
     2241            /*Vertex 1*/
     2242            xyz_list_RGB[3*1+0]=xyz_list[3*2+0]+(xyz_list[3*0+0]-xyz_list[3*2+0])*fraction1; //x
     2243            xyz_list_RGB[3*1+1]=xyz_list[3*2+1]+(xyz_list[3*0+1]-xyz_list[3*2+1])*fraction1; //y
     2244            /*Vertex 2*/
     2245            xyz_list_RGB[3*2+0]=xyz_list[3*2+0]+(xyz_list[3*1+0]-xyz_list[3*2+0])*fraction2; //x;
     2246            xyz_list_RGB[3*2+1]=xyz_list[3*2+1]+(xyz_list[3*1+1]-xyz_list[3*2+1])*fraction2; //y;
     2247            break;
     2248         default:
     2249            _error_("index "<<point1<<" not supported yet");
     2250      }//}}}
     2251   }
     2252        if(ig>3 && ig<8){ // GREEN element
     2253      switch(point1){//{{{
     2254         case 0:
     2255            /*Vertex 0 - local in the subelement*/
     2256            xyz_list_RGB[3*0+0]=xyz_list[3*0+0]+(xyz_list[3*1+0]-xyz_list[3*0+0])*fraction1; //x
     2257            xyz_list_RGB[3*0+1]=xyz_list[3*0+1]+(xyz_list[3*1+1]-xyz_list[3*0+1])*fraction1; //y
     2258            /*Vertex 1*/
     2259            xyz_list_RGB[3*1+0]=xyz_list[3*2+0]; //x
     2260            xyz_list_RGB[3*1+1]=xyz_list[3*2+1]; //y
     2261            /*Vertex 2*/
     2262            xyz_list_RGB[3*2+0]=xyz_list[3*0+0]+(xyz_list[3*2+0]-xyz_list[3*0+0])*fraction2; //x;
     2263            xyz_list_RGB[3*2+1]=xyz_list[3*0+1]+(xyz_list[3*2+1]-xyz_list[3*0+1])*fraction2; //y;
     2264            break;
     2265                         case 1:
     2266            /*Vertex 0*/
     2267            xyz_list_RGB[3*0+0]=xyz_list[3*1+0]+(xyz_list[3*2+0]-xyz_list[3*1+0])*fraction1; //x
     2268            xyz_list_RGB[3*0+1]=xyz_list[3*1+1]+(xyz_list[3*2+1]-xyz_list[3*1+1])*fraction1; //y
     2269            /*Vertex 1*/
     2270            xyz_list_RGB[3*1+0]=xyz_list[3*0+0]; //x
     2271            xyz_list_RGB[3*1+1]=xyz_list[3*0+1]; //y
     2272            /*Vertex 2*/
     2273            xyz_list_RGB[3*2+0]=xyz_list[3*1+0]+(xyz_list[3*0+0]-xyz_list[3*1+0])*fraction2; //x;
     2274            xyz_list_RGB[3*2+1]=xyz_list[3*1+1]+(xyz_list[3*0+1]-xyz_list[3*1+1])*fraction2; //y;
     2275            break;
     2276                        case 2:
     2277            /*Vertex 0*/
     2278            xyz_list_RGB[3*0+0]=xyz_list[3*2+0]+(xyz_list[3*0+0]-xyz_list[3*2+0])*fraction1; //x
     2279            xyz_list_RGB[3*0+1]=xyz_list[3*2+1]+(xyz_list[3*0+1]-xyz_list[3*2+1])*fraction1; //y
     2280            /*Vertex 1*/
     2281            xyz_list_RGB[3*1+0]=xyz_list[3*1+0]; //x
     2282            xyz_list_RGB[3*1+1]=xyz_list[3*1+1]; //y
     2283            /*Vertex 2*/
     2284            xyz_list_RGB[3*2+0]=xyz_list[3*2+0]+(xyz_list[3*1+0]-xyz_list[3*2+0])*fraction2; //x;
     2285            xyz_list_RGB[3*2+1]=xyz_list[3*2+1]+(xyz_list[3*1+1]-xyz_list[3*2+1])*fraction2; //y;
     2286            break;
     2287         default:
     2288            _error_("index "<<point1<<" not supported yet");
     2289      }//}}}
     2290   }
     2291        if(ig>7){ // RED element
     2292      switch(point1){//{{{
     2293         case 0:
     2294            /*Vertex 0*/
     2295            xyz_list_RGB[3*0+0]=xyz_list[3*0+0]+(xyz_list[3*1+0]-xyz_list[3*0+0])*fraction1; //x
     2296            xyz_list_RGB[3*0+1]=xyz_list[3*0+1]+(xyz_list[3*1+1]-xyz_list[3*0+1])*fraction1; //y
     2297            /*Vertex 1*/
     2298            xyz_list_RGB[3*1+0]=xyz_list[3*1+0]; //x
     2299            xyz_list_RGB[3*1+1]=xyz_list[3*1+1]; //y
     2300            /*Vertex 2*/
     2301            xyz_list_RGB[3*2+0]=xyz_list[3*2+0]; //x;
     2302            xyz_list_RGB[3*2+1]=xyz_list[3*2+1]; //y;
     2303            break;
     2304         case 1:
     2305                                /*Vertex 0*/
     2306            xyz_list_RGB[3*0+0]=xyz_list[3*1+0]+(xyz_list[3*2+0]-xyz_list[3*1+0])*fraction1; //x
     2307            xyz_list_RGB[3*0+1]=xyz_list[3*1+1]+(xyz_list[3*2+1]-xyz_list[3*1+1])*fraction1; //y
     2308            /*Vertex 1*/
     2309            xyz_list_RGB[3*1+0]=xyz_list[3*2+0]; //x
     2310            xyz_list_RGB[3*1+1]=xyz_list[3*2+1]; //y
     2311            /*Vertex 2*/
     2312            xyz_list_RGB[3*2+0]=xyz_list[3*0+0]; //x;
     2313            xyz_list_RGB[3*2+1]=xyz_list[3*0+1]; //y;
     2314            break;
     2315                        case 2:
     2316            /*Vertex 0*/
     2317            xyz_list_RGB[3*0+0]=xyz_list[3*2+0]+(xyz_list[3*0+0]-xyz_list[3*2+0])*fraction1; //x
     2318            xyz_list_RGB[3*0+1]=xyz_list[3*2+1]+(xyz_list[3*0+1]-xyz_list[3*2+1])*fraction1; //y
     2319            /*Vertex 1*/
     2320            xyz_list_RGB[3*1+0]=xyz_list[3*0+0]; //x
     2321            xyz_list_RGB[3*1+1]=xyz_list[3*0+1]; //y
     2322            /*Vertex 2*/
     2323            xyz_list_RGB[3*2+0]=xyz_list[3*1+0]; //x;
     2324            xyz_list_RGB[3*2+1]=xyz_list[3*1+1]; //y;
     2325            break;
     2326         default:
     2327            _error_("index "<<point1<<" not supported yet");
     2328      }//}}}
     2329   }
     2330
     2331        //std::cout<<"  ig="<<ig<<std::endl;
     2332        //std::cout<<"  "<<xyz_list_RGB[0*numnodes+0]<<"\t"<<xyz_list_RGB[1*numnodes+0]<<"\t"<<xyz_list_RGB[2*numnodes+0]<<std::endl; //x
     2333        //std::cout<<"  "<<xyz_list_RGB[0*numnodes+1]<<"\t"<<xyz_list_RGB[1*numnodes+1]<<"\t"<<xyz_list_RGB[2*numnodes+1]<<std::endl; //y
     2334
     2335   /*Get Jacobian*/
     2336        tria->GetJacobianInvert(&Jinv[0][0],xyz_list_RGB,gauss_CG); // gauss is not used here
     2337       
     2338        /*Build dbasis CG in the subelement*/
     2339        for(int i=0;i<numnodes;i++){
     2340                dbasis_subelem[numnodes*0+i]=Jinv[0][0]*dbasis_ref[0*numnodes+i]+Jinv[0][1]*dbasis_ref[1*numnodes+i];//dPhi_i/dx
     2341                dbasis_subelem[numnodes*1+i]=Jinv[1][0]*dbasis_ref[0*numnodes+i]+Jinv[1][1]*dbasis_ref[1*numnodes+i];//dPhi_i/dy
     2342        }
     2343
     2344   xDelete<IssmDouble>(xyz_list);
     2345   xDelete<IssmDouble>(xyz_list_RGB);
     2346   //xDelete<IssmDouble>(dbasis_subelem);
     2347
     2348}/*}}}*/
     2349#endif
    20352350
    20362351/*L1L2*/

  • issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h

    r25379 r25451  
    4545                void           GetBSSAprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
    4646                void           InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element);
     47                void                            ComputeSurfaceSlope(IssmDouble* slope,Gauss* gauss_DG,Gauss* gauss_CG,int point1,IssmDouble fraction1,IssmDouble fraction2,int ig,int dim,Element* element);
     48                void                            NodalFunctionsDerivativesRGB(IssmDouble* dbasis_subelem,Gauss* gauss_DG,Gauss* gauss_CG,int point1,IssmDouble fraction1,IssmDouble fraction2,int ig,int dim,Element* element);
     49
    4750                /*L1L2*/
    4851                ElementMatrix* CreateKMatrixL1L2(Element* element);

  • issm/trunk-jpl/src/c/classes/gauss/GaussTria.cpp

    r25446 r25451  
    366366         *
    367367         */
    368         int         ig;
     368        int         ii;
    369369        IssmDouble x,y;
    370370        IssmDouble xy_list[3][2];
    … …  
    384384        this->weights=xNew<IssmDouble>(this->numgauss);
    385385
    386         for(ig=0;ig<gauss1->numgauss;ig++){ // Add the first triangle gauss points (BLUE)
    387                 this->coords1[ig]=gauss1->coords1[ig];
    388                 this->coords2[ig]=gauss1->coords2[ig];
    389                 this->coords3[ig]=gauss1->coords3[ig];
    390                 this->weights[ig]=gauss1->weights[ig]*r1*r2;
    391         }
    392         for(ig=0;ig<gauss2->numgauss;ig++){ // Add the second triangle gauss points (GREEN)
    393                 this->coords1[gauss1->numgauss+ig]=gauss2->coords1[ig];
    394                 this->coords2[gauss1->numgauss+ig]=gauss2->coords2[ig];
    395                 this->coords3[gauss1->numgauss+ig]=gauss2->coords3[ig];
    396                 this->weights[gauss1->numgauss+ig]=gauss2->weights[ig]*r1*(1-r2);
    397         }
    398         for(ig=0;ig<gauss3->numgauss;ig++){ // Add the second triangle gauss points (RED)
    399                 this->coords1[gauss1->numgauss+gauss2->numgauss+ig]=gauss3->coords1[ig];
    400                 this->coords2[gauss1->numgauss+gauss2->numgauss+ig]=gauss3->coords2[ig];
    401                 this->coords3[gauss1->numgauss+gauss2->numgauss+ig]=gauss3->coords3[ig];
    402                 this->weights[gauss1->numgauss+gauss2->numgauss+ig]=gauss3->weights[ig]*(1-r1);
     386        for(ii=0;ii<gauss1->numgauss;ii++){ // Add the first triangle gauss points (BLUE)
     387                this->coords1[ii]=gauss1->coords1[ii];
     388                this->coords2[ii]=gauss1->coords2[ii];
     389                this->coords3[ii]=gauss1->coords3[ii];
     390                this->weights[ii]=gauss1->weights[ii]*r1*r2;
     391        }
     392        for(ii=0;ii<gauss2->numgauss;ii++){ // Add the second triangle gauss points (GREEN)
     393                this->coords1[gauss1->numgauss+ii]=gauss2->coords1[ii];
     394                this->coords2[gauss1->numgauss+ii]=gauss2->coords2[ii];
     395                this->coords3[gauss1->numgauss+ii]=gauss2->coords3[ii];
     396                this->weights[gauss1->numgauss+ii]=gauss2->weights[ii]*r1*(1-r2);
     397        }
     398        for(ii=0;ii<gauss3->numgauss;ii++){ // Add the second triangle gauss points (RED)
     399                this->coords1[gauss1->numgauss+gauss2->numgauss+ii]=gauss3->coords1[ii];
     400                this->coords2[gauss1->numgauss+gauss2->numgauss+ii]=gauss3->coords2[ii];
     401                this->coords3[gauss1->numgauss+gauss2->numgauss+ii]=gauss3->coords3[ii];
     402                this->weights[gauss1->numgauss+gauss2->numgauss+ii]=gauss3->weights[ii]*(1-r1);
    403403        }
    404404
Note: See TracChangeset for help on using the changeset viewer.