Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5650) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5651) @@ -2384,26 +2384,7 @@ /* 3d gaussian points: */ - int num_gauss,ig; - double* first_gauss_area_coord = NULL; - double* second_gauss_area_coord = NULL; - double* third_gauss_area_coord = NULL; - double* fourth_gauss_vert_coord = NULL; - double* area_gauss_weights = NULL; - double* vert_gauss_weights = NULL; - int ig1,ig2; - double gauss_weight1,gauss_weight2; - double gauss_coord[4]; - int order_area_gauss; - int num_vert_gauss; - int num_area_gauss; - double gauss_weight; - - /* 2d gaussian point: */ - int num_gauss2d; - double* first_gauss_area_coord2d = NULL; - double* second_gauss_area_coord2d = NULL; - double* third_gauss_area_coord2d = NULL; - double* gauss_weights2d=NULL; - double gauss_l1l2l3[3]; + int ig; + GaussPenta *gauss=NULL; + GaussTria *gauss_tria=NULL; /* material data: */ @@ -2471,15 +2452,4 @@ GetDofList(&doflistp,PattynApproximationEnum); //MacAyeal dof list - /*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.*/ - - order_area_gauss=5; - num_vert_gauss=5; - - gaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss); - /*Retrieve all inputs we will be needing: */ vx_input=inputs->GetInput(VxEnum); @@ -2489,55 +2459,42 @@ /* Start looping on the number of gaussian points: */ - for (ig1=0; ig1GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input); - this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss_coord,vxold_input,vyold_input); - - /*Get viscosity: */ - matice->GetViscosity3d(&viscosity, &epsilon[0]); - matice->GetViscosity3d(&oldviscosity, &oldepsilon[0]); - - /*Get B and Bprime matrices: */ - GetBMacAyealPattyn(&B[0][0], &xyz_list[0][0], gauss_coord); - tria->GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss_coord); - - /* Get Jacobian determinant: */ - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord); - - /*Build the D matrix: we plug the gaussian weight, the viscosity, and the jacobian determinant - onto this scalar matrix, so that we win some computational time: */ - - newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); - D_scalar=2*newviscosity*gauss_weight*Jdet; - for (i=0;i<3;i++){ - D[i][i]=D_scalar; - } - - /* Do the triple product tB*D*Bprime: */ - TripleMultiply( &B[0][0],3,numdofp,1, - &D[0][0],3,3,0, - &Bprime[0][0],3,numdofm,0, - &Ke_gg_gaussian[0][0],0); - - /* Add the Ke_gg_gaussian, and optionally Ke_gg_gaussian onto Ke_gg: */ - for( i=0; ibegin();igend();ig++){ + + gauss->GaussPoint(ig); + gauss->SynchronizeGaussTria(gauss_tria); + + /*Get strain rate from velocity: */ + this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input); + this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input); + + /*Get viscosity: */ + matice->GetViscosity3d(&viscosity, &epsilon[0]); + matice->GetViscosity3d(&oldviscosity, &oldepsilon[0]); + + /*Get B and Bprime matrices: */ + GetBMacAyealPattyn(&B[0][0], &xyz_list[0][0], gauss); + tria->GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss_tria); + + /* Get Jacobian determinant: */ + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + + /*Build the D matrix: we plug the gaussian weight, the viscosity, and the jacobian determinant + onto this scalar matrix, so that we win some computational time: */ + + newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); + D_scalar=2*newviscosity*gauss->weight*Jdet; + for (i=0;i<3;i++) D[i][i]=D_scalar; + + /* Do the triple product tB*D*Bprime: */ + TripleMultiply( &B[0][0],3,numdofp,1, + &D[0][0],3,3,0, + &Bprime[0][0],3,numdofm,0, + &Ke_gg_gaussian[0][0],0); + + /* Add the Ke_gg_gaussian, and optionally Ke_gg_gaussian onto Ke_gg: */ + for( i=0; imatice; delete tria; - xfree((void**)&first_gauss_area_coord); - xfree((void**)&second_gauss_area_coord); - xfree((void**)&third_gauss_area_coord); - xfree((void**)&fourth_gauss_vert_coord); - xfree((void**)&area_gauss_weights); - xfree((void**)&vert_gauss_weights); - xfree((void**)&first_gauss_area_coord2d); - xfree((void**)&second_gauss_area_coord2d); - xfree((void**)&third_gauss_area_coord2d); - xfree((void**)&gauss_weights2d); xfree((void**)&doflistm); xfree((void**)&doflistp); + delete gauss; + delete gauss_tria; } /*}}}*/ @@ -3113,18 +3062,6 @@ /* 3d gaussian points: */ - int num_gauss,ig; - double* first_gauss_area_coord = NULL; - double* second_gauss_area_coord = NULL; - double* third_gauss_area_coord = NULL; - double* fourth_gauss_vert_coord = NULL; - double* area_gauss_weights = NULL; - double* vert_gauss_weights = NULL; - int ig1,ig2; - double gauss_weight1,gauss_weight2; - double gauss_coord[4]; - int order_area_gauss; - int num_vert_gauss; - int num_area_gauss; - double gauss_weight; + int ig; + GaussPenta *gauss=NULL; /* matrices: */ @@ -3164,62 +3101,33 @@ GetDofList(&doflist); - /*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.*/ - - order_area_gauss=2; - num_vert_gauss=2; - - gaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss); - /* Start looping on the number of gaussian points: */ - for (ig1=0; ig1begin();igend();ig++){ + + gauss->GaussPoint(ig); + + /*Get B and Bprime matrices: */ + GetBVert(&B[0][0], &xyz_list[0][0], gauss); + GetBprimeVert(&Bprime[0][0], &xyz_list[0][0], gauss); + + /* Get Jacobian determinant: */ + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + DL_scalar=gauss->weight*Jdet; + + /* Do the triple product tB*D*Bprime: */ + TripleMultiply( &B[0][0],1,numgrids,1, + &DL_scalar,1,1,0, + &Bprime[0][0],1,numgrids,0, + &Ke_gg_gaussian[0][0],0); + + /* Add the Ke_gg_gaussian, and optionally Ke_gg_gaussian onto Ke_gg: */ + for( i=0; iparameters->FindParam(&epsvel,EpsVelEnum); - /* 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: */ - area_order=2; - num_vert_gauss=2; - - 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); + vx_input=inputs->GetInput(VxEnum); + vy_input=inputs->GetInput(VyEnum); + vz_input=inputs->GetInput(VzEnum); /* Start looping on the number of gaussian points: */ - for (igarea=0; igareaGetParameterValue(&u, gauss_coord,VxEnum); - inputs->GetParameterValue(&v, gauss_coord,VyEnum); - inputs->GetParameterValue(&w, gauss_coord,VzEnum); - - D_scalar=gauss_weight*Jdet; - - if(dt){ - D_scalar=D_scalar*dt; - } - - D[0][0]=D_scalar*u;D[0][1]=0; D[0][2]=0; - D[1][0]=0; D[1][1]=D_scalar*v;D[1][2]=0; - D[2][0]=0; D[2][1]=0; D[2][2]=D_scalar*w; - - /* Do the triple product B'*D*Bprime: */ - MatrixMultiply(&B_advec[0][0],3,numdof,1,&D[0][0],3,3,0,&tBD_advec[0][0],0); - MatrixMultiply(&tBD_advec[0][0],numdof,3,0,&Bprime_advec[0][0],3,numdof,0,&Ke_gaussian_advec[0][0],0); - - /*Transient: */ - if(dt){ - GetNodalFunctionsP1(&L[0], gauss_coord); - D_scalar=gauss_weight*Jdet; - D_scalar=D_scalar; - - /* Do the triple product L'*D*L: */ - MatrixMultiply(&L[0],numdof,1,0,&D_scalar,1,1,0,&tLD[0],0); - MatrixMultiply(&tLD[0],numdof,1,0,&L[0],1,numdof,0,&Ke_gaussian_transient[0][0],0); - } - else{ - for(i=0;ibegin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + + /*Conduction: */ + + /*Get B_conduct matrix: */ + GetBConduct(&B_conduct[0][0],&xyz_list[0][0],gauss); + + /*Build D: */ + D_scalar=gauss->weight*Jdet*(thermalconductivity/(rho_ice*heatcapacity)); + + if(dt) D_scalar=D_scalar*dt; + + D[0][0]=D_scalar; D[0][1]=0; D[0][2]=0; + D[1][0]=0; D[1][1]=D_scalar; D[1][2]=0; + D[2][0]=0; D[2][1]=0; D[2][2]=D_scalar; + + /* Do the triple product B'*D*B: */ + MatrixMultiply(&B_conduct[0][0],3,numdof,1,&D[0][0],3,3,0,&tBD_conduct[0][0],0); + MatrixMultiply(&tBD_conduct[0][0],numdof,3,0,&B_conduct[0][0],3,numdof,0,&Ke_gaussian_conduct[0][0],0); + + /*Advection: */ + + /*Get B_advec and Bprime_advec matrices: */ + GetBAdvec(&B_advec[0][0],&xyz_list[0][0],gauss); + GetBprimeAdvec(&Bprime_advec[0][0],&xyz_list[0][0],gauss); + + //Build the D matrix + vx_input->GetParameterValue(&u, gauss); + vy_input->GetParameterValue(&v, gauss); + vz_input->GetParameterValue(&w, gauss); + + D_scalar=gauss->weight*Jdet; + + if(dt) D_scalar=D_scalar*dt; + + D[0][0]=D_scalar*u;D[0][1]=0; D[0][2]=0; + D[1][0]=0; D[1][1]=D_scalar*v;D[1][2]=0; + D[2][0]=0; D[2][1]=0; D[2][2]=D_scalar*w; + + /* Do the triple product B'*D*Bprime: */ + MatrixMultiply(&B_advec[0][0],3,numdof,1,&D[0][0],3,3,0,&tBD_advec[0][0],0); + MatrixMultiply(&tBD_advec[0][0],numdof,3,0,&Bprime_advec[0][0],3,numdof,0,&Ke_gaussian_advec[0][0],0); + + /*Transient: */ + if(dt){ + GetNodalFunctionsP1(&L[0], gauss); + D_scalar=gauss->weight*Jdet; + D_scalar=D_scalar; + + /* Do the triple product L'*D*L: */ + MatrixMultiply(&L[0],numdof,1,0,&D_scalar,1,1,0,&tLD[0],0); + MatrixMultiply(&tLD[0],numdof,1,0,&L[0],1,numdof,0,&Ke_gaussian_transient[0][0],0); + } + else{ + for(i=0;iweight*Jdet/(pow(u,2)+pow(v,2)+epsvel); + if(dt) D_scalar=D_scalar*dt; + K[0][0]=D_scalar*pow(u,2); K[0][1]=D_scalar*fabs(u)*fabs(v); + K[1][0]=D_scalar*fabs(u)*fabs(v);K[1][1]=D_scalar*pow(v,2); + + /*Get B_artdiff: */ + GetBArtdiff(&B_artdiff[0][0],&xyz_list[0][0],gauss); + + /* Do the triple product B'*K*B: */ + MatrixMultiply(&B_artdiff[0][0],2,numdof,1,&K[0][0],2,2,0,&tBD_artdiff[0][0],0); + MatrixMultiply(&tBD_artdiff[0][0],numdof,2,0,&B_artdiff[0][0],2,numdof,0,&Ke_gaussian_artdiff[0][0],0); + } + else{ + for(i=0;iGetInput(ThicknessEnum); @@ -3959,55 +3802,40 @@ /* Start looping on the number of gaussian points: */ - for (ig1=0; ig1GetParameterValue(&thickness, gauss_coord); - - /*Compute slope at gaussian point: */ - surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss_coord); - - /* Get Jacobian determinant: */ - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord); - - /*Get nodal functions: */ - GetNodalFunctionsP1(l1l6, gauss_coord); - - /*Compute driving stress: */ - driving_stress_baseline=matpar->GetRhoIce()*matpar->GetG(); - - /*Build pe_g_gaussian vector: */ - for (i=0;ibegin();igend();ig++){ + + gauss->GaussPoint(ig); + + /*Compute thickness at gaussian point: */ + thickness_input->GetParameterValue(&thickness, gauss); + + /*Compute slope at gaussian point: */ + surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss); + + /* Get Jacobian determinant: */ + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + + /*Get nodal functions: */ + GetNodalFunctionsP1(l1l6, gauss); + + /*Compute driving stress: */ + driving_stress_baseline=matpar->GetRhoIce()*matpar->GetG(); + + /*Build pe_g_gaussian vector: */ + for (i=0;iweight*l1l6[i]; } - - /*Add pe_g_gaussian vector to pe_g: */ - for( i=0; iGetInput(VxEnum); @@ -4362,51 +4172,35 @@ /* Start looping on the number of gaussian points: */ - for (ig1=0; ig1GetParameterDerivativeValue(&du[0],&xyz_list[0][0],gauss_coord); - vy_input->GetParameterDerivativeValue(&dv[0],&xyz_list[0][0],gauss_coord); - dudx=du[0]; - dvdy=dv[1]; - - /* Get Jacobian determinant: */ - GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord); - - /*Get nodal functions: */ - GetNodalFunctionsP1(l1l6, gauss_coord); - - /*Build pe_g_gaussian vector: */ - for (i=0;ibegin();igend();ig++){ + + gauss->GaussPoint(ig); + + /*Get velocity derivative, with respect to x and y: */ + vx_input->GetParameterDerivativeValue(&du[0],&xyz_list[0][0],gauss); + vy_input->GetParameterDerivativeValue(&dv[0],&xyz_list[0][0],gauss); + dudx=du[0]; + dvdy=dv[1]; + + /* Get Jacobian determinant: */ + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + + /*Get nodal functions: */ + GetNodalFunctionsP1(l1l6, gauss); + + /*Build pe_g_gaussian vector: */ + for (i=0;iweight*l1l6[i]; + } + + /*Add pe_g_gaussian vector to pe_g: */ + for( i=0; iGetMeltingPoint(); - /* 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: */ - 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); vy_input=inputs->GetInput(VyEnum); vz_input=inputs->GetInput(VzEnum); + if (dt) temperature_input=inputs->GetInput(TemperatureEnum); /* 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 */ - GetNodalFunctionsP1(&L[0], gauss_coord); - - /*Build deformational heating: */ - GetPhi(&phi, &epsilon[0], viscosity); - - /*Build pe_gaussian */ - scalar_def=phi/(rho_ice*heatcapacity)*Jdet*gauss_weight; - if(dt){ - scalar_def=scalar_def*dt; - } - - for(i=0;iGetParameterValue(&temperature, gauss_coord,TemperatureEnum); - scalar_transient=temperature*Jdet*gauss_weight; - 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 */ + GetNodalFunctionsP1(&L[0], gauss); + + /*Build deformational heating: */ + GetPhi(&phi, &epsilon[0], viscosity); + + /*Build pe_gaussian */ + scalar_def=phi/(rho_ice*heatcapacity)*Jdet*gauss->weight; + if(dt) scalar_def=scalar_def*dt; + + for(i=0;iGetParameterValue(&temperature, gauss); + scalar_transient=temperature*Jdet*gauss->weight; + for(i=0;iGetVxStrainRate3d(epsilonvx,xyz_list,gauss); + vy_input->GetVyStrainRate3d(epsilonvy,xyz_list,gauss); + vz_input->GetVzStrainRate3d(epsilonvz,xyz_list,gauss); + + /*Sum all contributions*/ + for(i=0;i<6;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]+epsilonvz[i]; + +} +/*}}}*/ +/*FUNCTION Penta::GetStrainRate3d{{{1*/ +void Penta::GetStrainRate3d(double* epsilon,double* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input, Input* vz_input){ /*Compute the 3d Strain Rate (6 components): * Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 5650) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 5651) @@ -170,4 +170,5 @@ void GetStrainRate3dPattyn(double* epsilon,double* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input); void GetStrainRate3d(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input, Input* vz_input); + void GetStrainRate3d(double* epsilon,double* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input, Input* vz_input); Penta* GetUpperElement(void); Penta* GetLowerElement(void); Index: /issm/trunk/src/c/objects/Gauss/GaussPenta.cpp =================================================================== --- /issm/trunk/src/c/objects/Gauss/GaussPenta.cpp (revision 5650) +++ /issm/trunk/src/c/objects/Gauss/GaussPenta.cpp (revision 5651) @@ -225,2 +225,11 @@ } /*}}}*/ +/*FUNCTION GaussPenta::SynchronizeGaussTria{{{1*/ +void GaussPenta::SynchronizeGaussTria(GaussTria* gauss_tria){ + + gauss_tria->coord1=this->coord1; + gauss_tria->coord2=this->coord2; + gauss_tria->coord3=this->coord3; + gauss_tria->weight=UNDEF; +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Gauss/GaussPenta.h =================================================================== --- /issm/trunk/src/c/objects/Gauss/GaussPenta.h (revision 5650) +++ /issm/trunk/src/c/objects/Gauss/GaussPenta.h (revision 5651) @@ -9,4 +9,5 @@ /*{{{1*/ #include "./../../shared/shared.h" +class GaussTria; /*}}}*/ @@ -42,4 +43,5 @@ void GaussVertex(int iv); void GaussCenter(void); + void SynchronizeGaussTria(GaussTria* gauss_tria); }; #endif