Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5677) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5678) @@ -2807,5 +2807,4 @@ /* Start looping on the number of gaussian points: */ - /* Start looping on the number of gaussian points: */ gauss=new GaussPenta(5,5); for (ig=gauss->begin();igend();ig++){ @@ -3728,20 +3727,7 @@ /* gaussian points: */ - int num_area_gauss; - int igarea,igvert; - double* first_gauss_area_coord = NULL; - double* second_gauss_area_coord = NULL; - double* third_gauss_area_coord = NULL; - double* vert_gauss_coord = NULL; - double* area_gauss_weights = NULL; - double* vert_gauss_weights = NULL; - - /* specific gaussian point: */ - double gauss_weight,area_gauss_weight,vert_gauss_weight; - double gauss_coord[4]; - double gauss_coord_tria[3]; - - int area_order=5; - int num_vert_gauss=5; + int ig; + GaussPenta *gauss=NULL; + GaussTria *gauss_tria=NULL; double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ @@ -3805,13 +3791,4 @@ GetDofList(&doflist,StokesApproximationEnum); - /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore - get tria gaussian points as well as segment gaussian points. For tria gaussian - points, order of integration is 2, because we need to integrate the product tB*D*B' - which is a polynomial of degree 3 (see GaussLegendreTria for more details). For segment gaussian - points, same deal, which yields 3 gaussian points.*/ - - /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ - gaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss); - /*Retrieve all inputs we will be needing: */ vx_input=inputs->GetInput(VxEnum); @@ -3821,67 +3798,46 @@ /* Start looping on the number of gaussian points: */ - for (igarea=0; igareaGetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input,vz_input); - matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); - - /* Get Jacobian determinant: */ - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord); - - /* Get nodal functions */ - GetNodalFunctionsMINI(&l1l7[0], gauss_coord); - - /* Build gaussian vector */ - for(i=0;ibegin();igend();ig++){ + + gauss->GaussPoint(ig); + + /*Compute strain rate and viscosity: */ + this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); + + /* Get Jacobian determinant: */ + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + + /* Get nodal functions */ + GetNodalFunctionsMINI(&l1l7[0], gauss); + + /* Build gaussian vector */ + for(i=0;iweight*l1l7[i]; + } + + /*Add Pe_gaussian to terms in Pe_temp. Watch out for column orientation from matlab: */ + for(i=0;i<27;i++) Pe_temp[i]+=Pe_gaussian[i]; + + /*Get B and Bprime matrices: */ + GetBStokes(&B[0][0],&xyz_list[0][0],gauss); + GetBprimeStokes(&B_prime[0][0],&xyz_list[0][0], gauss); + + /*Get bubble part of Bprime */ + for(i=0;i<8;i++) for(j=0;j<3;j++) B_prime_bubble[i][j]=B_prime[i][j+24]; + + /* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant + * onto this scalar matrix, so that we win some computational time: */ + D_scalar=gauss->weight*Jdet; + for (i=0;i<6;i++) D[i][i]=D_scalar*2*viscosity; + for (i=6;i<8;i++) D[i][i]=-D_scalar*stokesreconditioning; + + /* Do the triple product tB*D*Bprime: */ + MatrixMultiply(&B[0][0],8,27,1,&D[0][0],8,8,0,&tBD[0][0],0); + MatrixMultiply(&tBD[0][0],27,8,0,&B_prime_bubble[0][0],8,3,0,&Ke_gaussian[0][0],0); + + /*Add Ke_gaussian and Ke_gaussian to terms in pKe. Watch out for column orientation from matlab: */ + for(i=0;i<27;i++) for(j=0;j<3;j++) Ke_temp[i][j]+=Ke_gaussian[i][j]; } @@ -3895,27 +3851,22 @@ } - xfree((void**)&first_gauss_area_coord); xfree((void**)&second_gauss_area_coord); xfree((void**)&third_gauss_area_coord); xfree((void**)&area_gauss_weights); - GaussLegendreTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, 2); - /* Start looping on the number of gauss 2d (nodes on the bedrock) */ - for (igarea=0; igareaGaussFaceTria(0,1,2,2); + gauss_tria=new GaussTria(); + for (ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + gauss->SynchronizeGaussTria(gauss_tria); /*Get the Jacobian determinant */ - tria->GetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0], gauss_coord_tria); + tria->GetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0], gauss_tria); /* Get bed at gaussian point */ - bed_input->GetParameterValue(&bed, gauss_coord); + bed_input->GetParameterValue(&bed, gauss); /*Get L matrix : */ - GetNodalFunctionsP1(l1l6, gauss_coord); + GetNodalFunctionsP1(l1l6, gauss); /*Get water_pressure at gaussian point */ @@ -3927,5 +3878,5 @@ for(i=0;iweight*Jdet2d*l1l6[i]*bed_normal[j]; } } @@ -3940,10 +3891,6 @@ /*Free ressources:*/ - xfree((void**)&first_gauss_area_coord); - xfree((void**)&second_gauss_area_coord); - xfree((void**)&third_gauss_area_coord); - xfree((void**)&area_gauss_weights); - xfree((void**)&vert_gauss_coord); - xfree((void**)&vert_gauss_weights); + delete gauss; + delete gauss_tria; xfree((void**)&doflist);