9 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
16 #include "../classes.h"
17 #include "../Inputs2/TriaInput2.h"
18 #include "../Inputs2/PentaInput2.h"
19 #include "../Inputs2/ControlInput2.h"
20 #include "../Inputs2/DatasetInput2.h"
21 #include "../Inputs2/TransientInput2.h"
22 #include "../../shared/shared.h"
24 #include "../../modules/GiaDeflectionCorex/GiaDeflectionCorex.h"
30 #define NUMVERTICES1D 2
99 for(i=0;i<nanalyses;i++){
113 else tria->
hnodes[i] = NULL;
134 unsigned int num_nodes = 3;
135 tria->
nodes = xNew<Node*>(num_nodes);
136 for(i=0;i<num_nodes;i++)
if(this->
nodes[i]) tria->
nodes[i]=this->
nodes[i];
else tria->
nodes[i] = NULL;
138 else tria->
nodes = NULL;
146 void Tria::Marshall(
char** pmarshalled_data,
int* pmarshalled_data_size,
int marshall_direction){
174 this->
AddInput2(input_enum,values,interpolation_enum);
177 _error_(
"not implemented yet");
192 int* temp = xNew<int>(3);
196 switch(interpolation_enum){
223 switch(interpolation_enum){
241 if(num_inputs<1)
_error_(
"Cannot create a DatasetInput of size <1");
243 if(N!=num_inputs)
_error_(
"Sizes are not consistent");
253 for(
int i=0;i<num_inputs;i++){
254 for(
int j=0;j<
NUMVERTICES;j++) nodeinputs[j]=array[vertexsids[j]*N+i];
261 bool already =
false;
277 partition[i]=reCast<int>(qmu_part[offsetsid[i]]);
278 mean=mean+1.0/
NUMVERTICES*vertex_response[offsetdof[i]];
285 if (partition[i]==partition[j]){
306 IssmDouble B,sigma_max,sigma_max_floating,sigma_max_grounded,n;
340 vel=sqrt(vx*vx+vy*vy)+1.e-14;
344 this->
StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
347 Matrix2x2Eigen(&lambda1,&lambda2,&ex,&ey,epsilon[0],epsilon[2],epsilon[1]);
348 _assert_(!xIsNan<IssmDouble>(lambda1));
349 _assert_(!xIsNan<IssmDouble>(lambda2));
352 lambda1 =
max(lambda1,0.);
353 lambda2 =
max(lambda2,0.);
356 epse_2 = 1./2. *(lambda1*lambda1 + lambda2*lambda2);
357 sigma_vm[iv] = sqrt(3.) * B * pow(epse_2,1./(2.*n));
366 sigma_max = sigma_max_floating;
368 sigma_max = sigma_max_grounded;
375 calvingratex[iv]=vx*sigma_vm[iv]/sigma_max;
376 calvingratey[iv]=vy*sigma_vm[iv]/sigma_max;
378 calvingrate[iv] =sqrt(calvingratex[iv]*calvingratex[iv] + calvingratey[iv]*calvingratey[iv]);
397 IssmDouble water_height, bed,Ho,thickness,surface;
399 IssmDouble B, strainparallel, straineffective,n;
401 int crevasse_opening_stress;
436 strainrateparallel_input->
GetInputValue(&strainparallel,gauss);
437 strainrateeffective_input->
GetInputValue(&straineffective,gauss);
447 vel=sqrt(vx*vx+vy*vy)+1.e-14;
449 s1=(s_xx+s_yy)/2.+sqrt(pow((s_xx-s_yy)/2.,2)+pow(s_xy,2));
450 s2=(s_xx+s_yy)/2.-sqrt(pow((s_xx-s_yy)/2.,2)+pow(s_xy,2));
451 if(fabs(s2)>fabs(s1)){stmp=s2; s2=s1; s1=stmp;}
453 Ho = thickness - (rho_seawater/rho_ice) * (-bed);
456 if(crevasse_opening_stress==0){
457 surface_crevasse[iv] = B * strainparallel * pow(straineffective, ((1 / n)-1)) / (rho_ice * constant_g);
458 basal_crevasse[iv] = (rho_ice/(rho_seawater-rho_ice)) * (B * strainparallel * pow(straineffective,((1/n)-1)) / (rho_ice*constant_g) - Ho);
460 else if(crevasse_opening_stress==1){
461 surface_crevasse[iv] = s1 / (rho_ice*constant_g);
462 basal_crevasse[iv] = (rho_ice/(rho_seawater-rho_ice))* (s1/ (rho_ice*constant_g)-Ho);
466 if (surface_crevasse[iv]<0.) {
467 surface_crevasse[iv]=0.;
470 if (basal_crevasse[iv]<0.) basal_crevasse[iv]=0.;
471 if (bed>0.) basal_crevasse[iv] = 0.;
478 surface_crevasse[iv] = surface_crevasse[iv] + (rho_freshwater/rho_ice)*water_height;
479 crevasse_depth[iv] = surface_crevasse[iv] + (rho_freshwater/rho_ice)*water_height + basal_crevasse[iv];
522 strainperpendicular_input->
GetInputValue(&strainperpendicular,gauss);
527 if(strainparallel>0. && strainperpendicular>0. && bed<=0.){
528 calvingrate[iv]=propcoeff*strainparallel*strainperpendicular;
533 calvingratex[iv]=calvingrate[iv]*vx/(sqrt(vel)+1.e-14);
534 calvingratey[iv]=calvingrate[iv]*vy/(sqrt(vel)+1.e-14);
555 int domaintype,index1,index2;
569 if(gl[0]==0.) gl[0]=gl[0]+epsilon;
570 if(gl[1]==0.) gl[1]=gl[1]+epsilon;
571 if(gl[2]==0.) gl[2]=gl[2]+epsilon;
582 s1=gl[2]/(gl[2]-gl[1]);
583 s2=gl[2]/(gl[2]-gl[0]);
587 xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
588 xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
589 xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
590 xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
591 xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
592 xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
594 else if(gl[1]*gl[2]>0){
597 s1=gl[0]/(gl[0]-gl[1]);
598 s2=gl[0]/(gl[0]-gl[2]);
603 xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
604 xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
605 xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
606 xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
607 xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
608 xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
610 else if(gl[0]*gl[2]>0){
613 s1=gl[1]/(gl[1]-gl[0]);
614 s2=gl[1]/(gl[1]-gl[2]);
619 xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
620 xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
621 xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
622 xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
623 xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
624 xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
640 normal[0] = -normal[0];
641 normal[1] = -normal[1];
646 IssmDouble calvingratex,calvingratey,thickness,Jdet,flux_per_area;
649 Input2* calvingratex_input=NULL;
650 Input2* calvingratey_input=NULL;
662 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
670 flux += rho_ice*Jdet*gauss->
weight*thickness*(calvingratex*normal[0] + calvingratey*normal[1]);
671 area += Jdet*gauss->
weight*thickness;
673 flux_per_area=flux/area;
691 int domaintype,index1,index2;
706 if(gl[0]==0.) gl[0]=gl[0]+epsilon;
707 if(gl[1]==0.) gl[1]=gl[1]+epsilon;
708 if(gl[2]==0.) gl[2]=gl[2]+epsilon;
719 s1=gl[2]/(gl[2]-gl[1]);
720 s2=gl[2]/(gl[2]-gl[0]);
724 xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
725 xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
726 xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
727 xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
728 xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
729 xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
731 else if(gl[1]*gl[2]>0){
734 s1=gl[0]/(gl[0]-gl[1]);
735 s2=gl[0]/(gl[0]-gl[2]);
740 xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
741 xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
742 xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
743 xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
744 xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
745 xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
747 else if(gl[0]*gl[2]>0){
750 s1=gl[1]/(gl[1]-gl[0]);
751 s2=gl[1]/(gl[1]-gl[2]);
756 xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
757 xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
758 xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
759 xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
760 xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
761 xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
777 normal[0] = -normal[0];
778 normal[1] = -normal[1];
783 IssmDouble calvingratex,calvingratey,vx,vy,vel,meltingrate,meltingratex,meltingratey,thickness,Jdet,flux_per_area;
786 Input2* calvingratex_input=NULL;
787 Input2* calvingratey_input=NULL;
790 Input2* meltingrate_input=NULL;
805 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
815 meltingratex=meltingrate*vx/(sqrt(vel)+1.e-14);
816 meltingratey=meltingrate*vy/(sqrt(vel)+1.e-14);
819 flux += rho_ice*Jdet*gauss->
weight*thickness*((calvingratex+meltingratex)*normal[0] + (calvingratey+meltingratey)*normal[1]);
820 area += Jdet*gauss->
weight*thickness;
822 flux_per_area=flux/area;
834 return sqrt(2*this->
GetArea());
838 _error_(
"Not Implemented yet");
856 int domaintype,dim=2;
876 this->
StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
877 switch(approximation){
892 tau_xx[iv]=2*viscosity*epsilon[0];
893 tau_yy[iv]=2*viscosity*epsilon[1];
894 tau_xy[iv]=2*viscosity*epsilon[2];
895 tau_e[iv]=1/sqrt(2)*sqrt(pow(tau_xx[iv],2)+pow(tau_yy[iv],2)+2*pow(tau_xy[iv],2));
898 Matrix2x2Eigen(&tau_2[iv],&tau_1[iv],NULL,NULL,tau_xx[iv],tau_xy[iv],tau_yy[iv]);
939 this->
StrainRateESA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
942 strain_xx[iv]=epsilon[0];
943 strain_yy[iv]=epsilon[1];
944 strain_xy[iv]=epsilon[2];
945 vorticity_xy[iv]=epsilon[3];
973 int domaintype,dim=2;
992 this->
StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
997 sigma_xx=2*viscosity*epsilon[0]-pressure;
998 sigma_yy=2*viscosity*epsilon[1]-pressure;
999 sigma_xy=2*viscosity*epsilon[2];
1005 sigma_nn[i]=sigma_xx*base_normal[0]*base_normal[0] + 2*sigma_xy*base_normal[0]*base_normal[1] + sigma_yy*base_normal[1]*base_normal[1];
1012 xDelete<IssmDouble>(xyz_list);
1013 xDelete<IssmDouble>(xyz_list_base);
1030 int domaintype,dim=2;
1048 this->
StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
1053 sigma_xx[iv]=2*viscosity*epsilon[0]-pressure;
1054 sigma_yy[iv]=2*viscosity*epsilon[1]-pressure;
1055 sigma_xy[iv]=2*viscosity*epsilon[2];
1073 int analysis_counter;
1083 else this->
hnodes[analysis_counter] = NULL;
1085 else this->
hnodes = NULL;
1091 else this->
nodes=NULL;
1113 for(
int n=0;n<N;n++){
1115 idlist[i] = offset + this->
vertices[i]->
Sid()+n*M;
1116 grad_list[i]=gradient[idlist[i]];
1130 for(
int i=0;i<
NUMVERTICES;i++) grad_list[i]=gradient[idlist[i]];
1161 values[iv] = reCast<IssmPDouble>(value)/reCast<IssmPDouble>(connectivity[iv]);
1162 gradients[iv] = reCast<IssmPDouble>(gradient)/reCast<IssmPDouble>(connectivity[iv]);
1179 if(xIsNan<IssmDouble>(distance))
_error_(
"NaN found in vector");
1180 if(xIsInf<IssmDouble>(distance))
_error_(
"Inf found in vector");
1199 int indices[3][2] = {{1,2},{2,0},{0,1}};
1204 for(
int i=0;i<3;i++){
1205 if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1.){
1210 _printf_(
"list of vertices on bed: "<<values[0]<<
" "<<values[1]<<
" "<<values[2]);
1211 _error_(
"Could not find 2 vertices on bed");
1217 int indices[3][2] = {{1,2},{2,0},{0,1}};
1222 for(
int i=0;i<3;i++){
1223 if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1.){
1224 *pindex1 = indices[i][0];
1225 *pindex2 = indices[i][1];
1230 _printf_(
"list of vertices on bed: "<<values[0]<<
" "<<values[1]<<
" "<<values[2]);
1231 _error_(
"Could not find 2 vertices on bed");
1237 int indices[3][2] = {{1,2},{2,0},{0,1}};
1242 for(
int i=0;i<3;i++){
1243 if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1.){
1248 _printf_(
"list of vertices on surface: "<<values[0]<<
" "<<values[1]<<
" "<<values[2]);
1249 _error_(
"Could not find 2 vertices on surface");
1255 int indices[3][2] = {{1,2},{2,0},{0,1}};
1260 for(
int i=0;i<3;i++){
1261 if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1.){
1262 *pindex1 = indices[i][0];
1263 *pindex2 = indices[i][1];
1268 _printf_(
"list of vertices on surface: "<<values[0]<<
" "<<values[1]<<
" "<<values[2]);
1269 _error_(
"Could not find 2 vertices on surface");
1274 switch(response_enum){
1303 xmin=xyz_list[0][0]; xmax=xyz_list[0][0];
1304 ymin=xyz_list[0][1]; ymax=xyz_list[0][1];
1307 if(xyz_list[i][0]<xmin) xmin=xyz_list[i][0];
1308 if(xyz_list[i][0]>xmax) xmax=xyz_list[i][0];
1309 if(xyz_list[i][1]<ymin) ymin=xyz_list[i][1];
1310 if(xyz_list[i][1]>ymax) ymax=xyz_list[i][1];
1337 floatingarea=(1-phi)*this->
GetArea();
1341 floatingarea=floatingarea*scalefactor;
1345 xDelete<IssmDouble>(xyz_list);
1346 return floatingarea;
1357 _error_(
" contact contiditon only works for FS elements");
1386 this->
StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
1393 sigmaxx[iv]=2*viscosity*epsilon[0]-pressure[iv];
1394 sigmayy[iv]=2*viscosity*epsilon[1]-pressure[iv];
1395 sigmaxy[iv]=2*viscosity*epsilon[2];
1403 sigma_nn[i]=-1*(sigmaxx[i]*bed_normal[0]*bed_normal[0] + sigmayy[i]*bed_normal[1]*bed_normal[1]+2*sigmaxy[i]*bed_normal[0]*bed_normal[1]);
1404 water_pressure[i]=-gravity*rho_water*base[i];
1420 xDelete<IssmDouble>(xyz_list);
1430 x1=xyz_list[0][0]; y1=xyz_list[0][1];
1431 x2=xyz_list[1][0]; y2=xyz_list[1][1];
1432 x3=xyz_list[2][0]; y3=xyz_list[2][1];
1434 _assert_(x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1>0);
1435 return (x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1)/2;
1451 x1=xyz_list[0][0]; y1=xyz_list[0][1]; z1=xyz_list[0][2];
1452 x2=xyz_list[1][0]; y2=xyz_list[1][1]; z2=xyz_list[1][2];
1453 x3=xyz_list[2][0]; y3=xyz_list[2][1]; z3=xyz_list[2][2];
1455 detm1=x1*y2 - x2*y1 - x1*y3 + x3*y1 + x2*y3 - x3*y2;
1456 detm2=y1*z2 - y2*z1 - y1*z3 + y3*z1 + y2*z3 - y3*z2;
1457 detm3=x2*z1 - x1*z2 + x1*z3 - x3*z1 - x2*z3 + x3*z2;
1459 return sqrt(pow(detm1,2) + pow(detm2,2) + pow(detm3,2))/2;
1473 switch (numiceverts){
1478 area_fraction=s[0]*s[1];
1481 area_fraction=s[0]+s[1]-s[0]*s[1];
1487 _error_(
"Wrong number of ice vertices in Tria::GetAreaIce!");
1490 _assert_((area_fraction>=0.) && (area_fraction<=1.));
1492 xDelete<int>(indices);
1493 return area_fraction*this->
GetArea();
1497 bool spherical=
true;
1500 IssmDouble arc12,arc23,arc31,semi_peri,excess;
1504 x1=llr_list[0][0]/180*
PI; y1=llr_list[0][1]/180*
PI; z1=llr_list[0][2];
1505 x2=llr_list[1][0]/180*
PI; y2=llr_list[1][1]/180*
PI; z2=llr_list[1][2];
1506 x3=llr_list[2][0]/180*
PI; y3=llr_list[2][1]/180*
PI; z3=llr_list[2][2];
1509 arc12=2.*asin(sqrt(pow(sin((x2-x1)/2),2.0)+cos(x1)*cos(x2)*pow(sin((y2-y1)/2),2)));
1510 arc23=2.*asin(sqrt(pow(sin((x3-x2)/2),2.0)+cos(x2)*cos(x3)*pow(sin((y3-y2)/2),2)));
1511 arc31=2.*asin(sqrt(pow(sin((x1-x3)/2),2.0)+cos(x3)*cos(x1)*pow(sin((y1-y3)/2),2)));
1514 semi_peri=(arc12+arc23+arc31)/2;
1517 excess=4.*atan(sqrt(tan(semi_peri/2)*tan((semi_peri-arc12)/2)*tan((semi_peri-arc23)/2)*tan((semi_peri-arc31)/2)));
1520 return excess*pow((z1+z2+z3)/3,2);
1535 xyz_bis[j][k]=xyz_list[j*3+k];
1538 for(i=0;i<numpoints;i++){
1542 xyz_bis[j][k]=xyz_zero[i*3+k];
1546 area_portion=fabs(xyz_bis[1][0]*xyz_bis[2][1] - xyz_bis[1][1]*xyz_bis[2][0] + xyz_bis[0][0]*xyz_bis[1][1] - xyz_bis[0][1]*xyz_bis[1][0] + xyz_bis[2][0]*xyz_bis[0][1] - xyz_bis[2][1]*xyz_bis[0][0])/2.;
1547 *(area_coordinates+3*i+j)=area_portion/area_init;
1552 xyz_bis[j][k]=xyz_list[j*3+k];
1578 if(gl[0]==0.) gl[0]=gl[0]+epsilon;
1579 if(gl[1]==0.) gl[1]=gl[1]+epsilon;
1580 if(gl[2]==0.) gl[2]=gl[2]+epsilon;
1583 if(gl[0]>0 && gl[1]>0 && gl[2]>0){
1588 else if(gl[0]<0 && gl[1]<0 && gl[2]<0){
1594 if(gl[0]*gl[1]*gl[2]<0) floating=
false;
1598 f1=gl[2]/(gl[2]-gl[0]);
1599 f2=gl[2]/(gl[2]-gl[1]);
1601 else if(gl[1]*gl[2]>0){
1603 f1=gl[0]/(gl[0]-gl[1]);
1604 f2=gl[0]/(gl[0]-gl[2]);
1606 else if(gl[0]*gl[2]>0){
1608 f1=gl[1]/(gl[1]-gl[2]);
1609 f2=gl[1]/(gl[1]-gl[0]);
1611 else _error_(
"case not possible");
1616 *mainlyfloating=floating;
1622 bool mainlyfloating =
true;
1623 int domaintype,index1,index2;
1625 IssmDouble phi,s1,s2,area_init,area_grounded;
1634 if(gl[0]==0.) gl[0]=gl[0]+epsilon;
1635 if(gl[1]==0.) gl[1]=gl[1]+epsilon;
1636 if(gl[2]==0.) gl[2]=gl[2]+epsilon;
1640 if(gl[index1]>0 && gl[index2]>0) phi=1;
1641 else if(gl[index1]<0 && gl[index2]<0) phi=0;
1642 else if(gl[index1]<0 && gl[index2]>0){
1643 phi=1./(1.-gl[index1]/gl[index2]);
1645 else if(gl[index2]<0 && gl[index1]>0){
1646 phi=1./(1.-gl[index2]/gl[index1]);
1652 if(gl[0]>0 && gl[1]>0 && gl[2]>0){
1655 else if(gl[0]<0 && gl[1]<0 && gl[2]<0){
1660 if(gl[0]*gl[1]*gl[2]>0) mainlyfloating=
false;
1664 xyz_bis[2][0]=xyz_list[3*2+0];
1665 xyz_bis[2][1]=xyz_list[3*2+1];
1666 xyz_bis[2][2]=xyz_list[3*2+2];
1669 s1=gl[2]/(gl[2]-gl[1]);
1670 s2=gl[2]/(gl[2]-gl[0]);
1673 xyz_bis[1][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
1674 xyz_bis[1][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
1675 xyz_bis[1][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
1678 xyz_bis[0][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
1679 xyz_bis[0][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
1680 xyz_bis[0][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
1682 else if(gl[1]*gl[2]>0){
1684 xyz_bis[0][0]=xyz_list[3*0+0];
1685 xyz_bis[0][1]=xyz_list[3*0+1];
1686 xyz_bis[0][2]=xyz_list[3*0+2];
1689 s1=gl[0]/(gl[0]-gl[1]);
1690 s2=gl[0]/(gl[0]-gl[2]);
1693 xyz_bis[1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
1694 xyz_bis[1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
1695 xyz_bis[1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
1698 xyz_bis[2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
1699 xyz_bis[2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
1700 xyz_bis[2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
1702 else if(gl[0]*gl[2]>0){
1704 xyz_bis[1][0]=xyz_list[3*1+0];
1705 xyz_bis[1][1]=xyz_list[3*1+1];
1706 xyz_bis[1][2]=xyz_list[3*1+2];
1709 s1=gl[1]/(gl[1]-gl[0]);
1710 s2=gl[1]/(gl[1]-gl[2]);
1713 xyz_bis[0][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
1714 xyz_bis[0][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
1715 xyz_bis[0][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
1718 xyz_bis[2][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
1719 xyz_bis[2][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
1720 xyz_bis[2][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
1722 else _error_(
"case not possible");
1733 if(mainlyfloating==
true) area_grounded=area_init-area_grounded;
1734 phi=area_grounded/area_init;
1739 if(phi>1 || phi<0)
_error_(
"Error. Problem with portion of grounded element: value should be between 0 and 1");
1752 int nrfrontbed,numiceverts;
1765 if(bed[i]<0.) nrfrontbed++;
1780 for(
int i=0;i<numiceverts;i++){
1782 x[counter] = xyz_list[indices[i]][0]+s[counter]*(xyz_list[indices[n]][0]-xyz_list[indices[i]][0]);
1783 y[counter] = xyz_list[indices[i]][1]+s[counter]*(xyz_list[indices[n]][1]-xyz_list[indices[i]][1]);
1791 if(lsf[indices[i]]==0.){
1792 x[counter]=xyz_list[indices[i]][0];
1793 y[counter]=xyz_list[indices[i]][1];
1796 if(counter==2)
break;
1806 _error_(
"not sure what's going on here...");
1808 x1=x[0]; y1=y[0]; x2=x[1]; y2=y[1];
1809 distance=sqrt(pow((x1-x2),2)+pow((y1-y2),2));
1811 int numthk=numiceverts+2;
1812 H=xNew<IssmDouble>(numthk);
1813 for(
int iv=0;iv<
NUMVERTICES;iv++) Haux[iv]=-bed[indices[iv]];
1815 switch(numiceverts){
1818 H[1]=Haux[0]+s[0]*(Haux[1]-Haux[0]);
1819 H[2]=Haux[0]+s[1]*(Haux[2]-Haux[0]);
1820 Haverage=(H[1]+H[2])/2;
1825 H[2]=Haux[0]+s[0]*(Haux[2]-Haux[0]);
1826 H[3]=Haux[1]+s[1]*(Haux[2]-Haux[1]);
1827 Haverage=(H[2]+H[3])/2;
1830 _error_(
"Number of ice covered vertices wrong in Tria::GetIceFrontArea(void)");
1833 frontarea=distance*Haverage;
1837 xDelete<int>(indices);
1838 xDelete<IssmDouble>(H);
1854 int num_frontnodes=0;
1856 if(levelset[i]>=0.){
1857 indicesfront[num_frontnodes]=i;
1865 int index=indicesfront[0];
1866 indicesfront[0]=indicesfront[1];
1867 indicesfront[1]=index;
1870 IssmDouble* xyz_front = xNew<IssmDouble>(3*2);
1872 for(
int i=0;i<2;i++){
1873 for(
int j=0;j<3;j++){
1874 xyz_front[3*i+j]=xyz_list[3*indicesfront[i]+j];
1878 *pxyz_front=xyz_front;
1885 if(!input)
return input;
1892 if(!input)
return input;
1903 if(!input)
return input;
1910 if(!input)
return input;
1921 if(!input)
return input;
1928 if(!input)
return input;
1948 for(
int iv=0;iv<
NUMVERTICES;iv++) pvalue[iv] = default_value;
1964 for(
int iv=0;iv<numnodes;iv++){
1970 for(
int iv=0;iv<numnodes;iv++) pvalue[iv] = default_value;
1977 if(!input_in)
return;
1982 PentaInput2* input = xDynamicCast<PentaInput2*>(input_in);
1990 switch(interpolation){
1993 indices[0] = this->
lid;
1994 input->
Serve(numindices,&indices[0]);
1998 for(
int i=0;i<3;i++) indices[i] =
vertices[i]->
lid;
1999 input->
Serve(numindices,&indices[0]);
2014 TriaInput2* input = xDynamicCast<TriaInput2*>(input_in);
2022 switch(interpolation){
2025 indices[0] = this->
lid;
2026 input->
Serve(numindices,&indices[0]);
2030 for(
int i=0;i<3;i++) indices[i] =
vertices[i]->
lid;
2031 input->
Serve(numindices,&indices[0]);
2035 input->
Serve(this->
lid,numindices);
2047 if(!datasetinput)
return NULL;
2049 for(
int i=0;i<datasetinput->
GetNumIds();i++){
2062 switch(interpolation){
2065 indices[0] = this->
lid;
2066 input->
Serve(numindices,&indices[0]);
2070 for(
int i=0;i<3;i++) indices[i] =
vertices[i]->
lid;
2071 input->
Serve(numindices,&indices[0]);
2093 switch(interpolation){
2096 indices[0] = this->
lid;
2097 input->
Serve(numindices,&indices[0]);
2101 for(
int i=0;i<3;i++) indices[i] =
vertices[i]->
lid;
2102 input->
Serve(numindices,&indices[0]);
2106 input->
Serve(this->
lid,numindices);
2114 return datasetinput;
2123 Input2* averaged_copy = averaged_input->
copy();
2125 averaged_copy->
ChangeEnum(averagedinput_enum);
2136 transient_input->
GetAllTimes(×teps,&numtimesteps);
2140 bool iscurrenttime_included =
false;
2141 for(
int i=0;i<numtimesteps;i++){
2142 if(timesteps[i]==currenttime) iscurrenttime_included=
true;
2143 if(timesteps[i]>currenttime)
break;
2146 if(iscurrenttime_included==
false)numsteps++;
2149 IssmDouble* times=xNew<IssmDouble>(numsteps);
2150 IssmDouble* values=xNew<IssmDouble>(numsteps);
2152 for(
int i=0;i<numsteps;i++){
2153 if((iscurrenttime_included==
false) && (i==(numsteps-1))){
2156 times[i]=currenttime;
2162 times[i]=timesteps[i];
2169 *pnumtimes=numtimesteps;
2203 int i, dir,nrfrontnodes;
2213 if(levelset[i]==level){
2214 indicesfront[nrfrontnodes]=i;
2223 int index=indicesfront[0];
2224 indicesfront[0]=indicesfront[1];
2225 indicesfront[1]=index;
2228 IssmDouble* xyz_front = xNew<IssmDouble>(3*nrfrontnodes);
2230 for(i=0;i<nrfrontnodes;i++){
2231 for(dir=0;dir<3;dir++){
2232 xyz_front[3*i+dir]=xyz_list[3*indicesfront[i]+dir];
2236 *pxyz_front=xyz_front;
2246 int i, numiceverts, numnoiceverts;
2247 int ind0, ind1, lastindex;
2261 if((lsf[ind0]-level)*(lsf[ind1]-level)<=0.){
2262 if(lsf[ind1]==level)
2275 indices_ice[numiceverts]=i;
2279 indices_noice[numnoiceverts]=i;
2284 for(i=0;i<numiceverts;i++){indices[i]=indices_ice[i];}
2285 for(i=0;i<numnoiceverts;i++){indices[numiceverts+i]=indices_noice[i];}
2287 switch (numiceverts){
2294 fraction[i]=(level-lsf[indices[0]])/(lsf[indices[numiceverts+i]]-lsf[indices[0]]);
2299 fraction[i]=(level-lsf[indices[i]])/(lsf[indices[numiceverts]]-lsf[indices[i]]);
2307 _error_(
"Wrong number of ice vertices in Tria::GetLevelsetIntersection!");
2312 *pnumiceverts=numiceverts;
2325 if(gl[0]==0.) gl[0]=gl[0]+epsilon;
2326 if(gl[1]==0.) gl[1]=gl[1]+epsilon;
2327 if(gl[2]==0.) gl[2]=gl[2]+epsilon;
2330 if(gl[0]>0 && gl[1]>0 && gl[2]>0){
2335 else if(gl[0]<0 && gl[1]<0 && gl[2]<0){
2341 if(gl[0]*gl[1]*gl[2]<0) negative=
false;
2345 f1=gl[2]/(gl[2]-gl[0]);
2346 f2=gl[2]/(gl[2]-gl[1]);
2348 else if(gl[1]*gl[2]>0){
2350 f1=gl[0]/(gl[0]-gl[1]);
2351 f2=gl[0]/(gl[0]-gl[2]);
2353 else if(gl[0]*gl[2]>0){
2355 f1=gl[1]/(gl[1]-gl[2]);
2356 f2=gl[1]/(gl[1]-gl[0]);
2359 _error_(
"This case should NOT be happening");
2365 *mainlynegative=negative;
2374 _error_(
"Vertex provided not found among element vertices");
2408 switch(interpolation){
2412 input->
Serve(numindices,&indices[0]);
2419 xDynamicCast<PentaInput2*>(input)->SetServeCollapsed(
true);
2450 TriaInput2* triainput = xDynamicCast<TriaInput2*>(input);
2463 TransientInput2* transientinput = xDynamicCast<TransientInput2*>(input);
2483 xDelete<int>(idlist);
2484 xDelete<IssmDouble>(values);
2500 IssmDouble* xyz_list_edge = xNew<IssmDouble>(2*3);
2502 for(
int i=0;i<2;i++)
for(
int j=0;j<2;j++) xyz_list_edge[i*3+j]=xyz_list[indices[i]][j];
2505 *pxyz_list = xyz_list_edge;
2517 IssmDouble* xyz_list_edge = xNew<IssmDouble>(2*3);
2519 for(
int i=0;i<2;i++)
for(
int j=0;j<2;j++) xyz_list_edge[i*3+j]=xyz_list[indices[i]][j];
2522 *pxyz_list = xyz_list_edge;
2540 groundedarea=phi*this->
GetArea();
2544 groundedarea=groundedarea*scalefactor;
2548 xDelete<IssmDouble>(xyz_list);
2549 return groundedarea;
2559 sum = values[0]+values[1]+values[2];
2561 _assert_(sum==0. || sum==1. || sum==2.);
2563 if(sum==3.)
_error_(
"Two edges on bed not supported yet...");
2580 sum = values[0]+values[1]+values[2];
2582 _assert_(sum==0. || sum==1. || sum==2.);
2584 if(sum==3.)
_error_(
"Two edges on surface not supported yet...");
2636 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
2644 flux += rho_ice*Jdet*gauss->
weight*thickness*(vx*normal[0] + vy*normal[1]);
2658 int domaintype,index1,index2;
2673 if(gl[0]==0.) gl[0]=gl[0]+epsilon;
2674 if(gl[1]==0.) gl[1]=gl[1]+epsilon;
2675 if(gl[2]==0.) gl[2]=gl[2]+epsilon;
2686 s1=gl[2]/(gl[2]-gl[1]);
2687 s2=gl[2]/(gl[2]-gl[0]);
2691 xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
2692 xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
2693 xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
2694 xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
2695 xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
2696 xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
2698 else if(gl[1]*gl[2]>0){
2701 s1=gl[0]/(gl[0]-gl[1]);
2702 s2=gl[0]/(gl[0]-gl[2]);
2707 xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
2708 xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
2709 xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
2710 xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
2711 xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
2712 xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
2714 else if(gl[0]*gl[2]>0){
2717 s1=gl[1]/(gl[1]-gl[0]);
2718 s2=gl[1]/(gl[1]-gl[2]);
2723 xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
2724 xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
2725 xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
2726 xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
2727 xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
2728 xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
2744 normal[0] = -normal[0];
2745 normal[1] = -normal[1];
2765 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
2773 flux += rho_ice*Jdet*gauss->
weight*thickness*(vx*normal[0] + vy*normal[1]);
2788 int domaintype,index1,index2;
2802 if(gl[0]==0.) gl[0]=gl[0]+epsilon;
2803 if(gl[1]==0.) gl[1]=gl[1]+epsilon;
2804 if(gl[2]==0.) gl[2]=gl[2]+epsilon;
2815 s1=gl[2]/(gl[2]-gl[1]);
2816 s2=gl[2]/(gl[2]-gl[0]);
2820 xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
2821 xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
2822 xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
2823 xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
2824 xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
2825 xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
2827 else if(gl[1]*gl[2]>0){
2830 s1=gl[0]/(gl[0]-gl[1]);
2831 s2=gl[0]/(gl[0]-gl[2]);
2836 xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
2837 xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
2838 xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
2839 xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
2840 xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
2841 xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
2843 else if(gl[0]*gl[2]>0){
2846 s1=gl[1]/(gl[1]-gl[0]);
2847 s2=gl[1]/(gl[1]-gl[2]);
2852 xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
2853 xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
2854 xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
2855 xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
2856 xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
2857 xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
2892 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
2900 flux += rho_ice*Jdet*gauss->
weight*thickness*(vx*normal[0] + vy*normal[1]);
2912 IssmDouble area_base,surface,base,Haverage,scalefactor;
2929 int numthk=numiceverts+2;
2930 H=xNew<IssmDouble>(numthk);
2935 for(i=0;i<
NUMVERTICES;i++) SFaux[i]= scalefactors[indices[i]];
2936 switch(numiceverts){
2939 SF[1]=SFaux[0]+s[0]*(SFaux[1]-SFaux[0]);
2940 SF[2]=SFaux[0]+s[1]*(SFaux[2]-SFaux[0]);
2945 SF[2]=SFaux[0]+s[0]*(SFaux[2]-SFaux[0]);
2946 SF[3]=SFaux[1]+s[1]*(SFaux[2]-SFaux[1]);
2949 _error_(
"Number of ice covered vertices wrong in Tria::IceVolume()");
2953 for(i=0;i<numthk;i++) scalefactor+=SF[i];
2955 area_base=area_base*scalefactor;
2959 for(i=0;i<
NUMVERTICES;i++) Haux[i]= surfaces[indices[i]]-bases[indices[i]];
2960 switch(numiceverts){
2963 H[1]=Haux[0]+s[0]*(Haux[1]-Haux[0]);
2964 H[2]=Haux[0]+s[1]*(Haux[2]-Haux[0]);
2969 H[2]=Haux[0]+s[0]*(Haux[2]-Haux[0]);
2970 H[3]=Haux[1]+s[1]*(Haux[2]-Haux[1]);
2973 _error_(
"Number of ice covered vertices wrong in Tria::IceVolume()");
2977 for(i=0;i<numthk;i++) Haverage+=H[i];
2986 area_base=area_base*scalefactor;
2994 Haverage=surface-base;
2998 xDelete<int>(indices);
2999 xDelete<IssmDouble>(H);
3000 xDelete<IssmDouble>(SF);
3006 return area_base*Haverage;
3014 IssmDouble base,surface,bed,bathymetry,scalefactor;
3026 base = 1./2. * fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
3030 base=base*scalefactor;
3037 if(!bed_input)
_error_(
"Could not find bed");
3043 return base*(surface-bed+
min(rho_water/rho_ice*bathymetry,0.));
3055 newinput=oldinput->
copy();
3068 int tria_vertex_ids[3];
3072 bool control_analysis,ad_analysis =
false;
3073 int num_control_type,num_responses;
3074 char** controls = NULL;
3079 iomodel->
FindConstant(&control_analysis,
"md.inversion.iscontrol");
3080 iomodel->
FindConstant(&ad_analysis,
"md.autodiff.isautodiff");
3081 if(control_analysis && !ad_analysis) iomodel->
FindConstant(&num_control_type,
"md.inversion.num_control_parameters");
3082 if(control_analysis && !ad_analysis) iomodel->
FindConstant(&num_responses,
"md.inversion.num_cost_functions");
3083 if(control_analysis && ad_analysis) iomodel->
FindConstant(&num_control_type,
"md.autodiff.num_independent_objects");
3084 if(control_analysis && ad_analysis) iomodel->
FindConstant(&num_responses,
"md.autodiff.num_dependent_objects");
3088 tria_vertex_ids[i]=reCast<int>(iomodel->
elements[3*index+i]);
3095 int* doflist = NULL;
3102 IssmDouble* values = xNew<IssmDouble>(numnodes);
3105 for(
int i=0;i<numnodes;i++){
3106 values[i]=solution[doflist[i]];
3107 if(xIsNan<IssmDouble>(values[i]))
_error_(
"NaN found in solution vector");
3108 if(xIsInf<IssmDouble>(values[i]))
_error_(
"Inf found in solution vector, SID = " << this->
Sid());
3115 xDelete<IssmDouble>(values);
3116 xDelete<int>(doflist);
3127 int *doflist = NULL;
3137 if(xIsNan<IssmDouble>(values[i]))
_error_(
"NaN found in vector");
3138 if(xIsInf<IssmDouble>(values[i]))
_error_(
"Inf found in vector");
3148 if(xIsNan<IssmDouble>(values[i]))
_error_(
"NaN found in vector");
3149 if(xIsInf<IssmDouble>(values[i]))
_error_(
"Inf found in vector");
3159 if(xIsNan<IssmDouble>(values[i]))
_error_(
"NaN found in vector");
3160 if(xIsInf<IssmDouble>(values[i]))
_error_(
"Inf found in vector");
3169 values = xNew<IssmDouble>(numnodes);
3172 for(
int i=0;i<numnodes;i++){
3173 values[i]=vector[doflist[i]];
3174 if(xIsNan<IssmDouble>(values[i]))
_error_(
"NaN found in vector");
3175 if(xIsInf<IssmDouble>(values[i]))
_error_(
"Inf found in vector");
3184 values = xNew<IssmDouble>(numnodes);
3186 for(
int i=0;i<numnodes;i++){
3188 if(xIsNan<IssmDouble>(values[i]))
_error_(
"NaN found in vector");
3189 if(xIsInf<IssmDouble>(values[i]))
_error_(
"Inf found in vector");
3200 value=vector[this->
Sid()];
3201 if(xIsNan<IssmDouble>(value))
_error_(
"NaN found in vector");
3202 if(xIsInf<IssmDouble>(value))
_error_(
"Inf found in vector");
3214 xDelete<int>(doflist);
3215 xDelete<IssmDouble>(values);
3226 sum = values[0]+values[1]+values[2];
3228 _assert_(sum==0. || sum==1. || sum==2.);
3230 if(sum==3.)
_error_(
"Two edges on boundary not supported yet...");
3253 if(ls[i]<0.) nrice++;
3254 if(nrice==1) isicefront=
true;
3283 if(ls[0]*ls[1]<0. || ls[0]*ls[2]<0. || (ls[0]*ls[1]*ls[2]==0. && ls[0]*ls[1]+ls[0]*ls[2]+ls[1]*ls[2]<=0.)){
3323 Input2* thickness_input=NULL;
3332 bool mainlynegative=
true;
3354 Gauss* gauss = this->
NewGauss(point1,fraction1,fraction2,mainlynegative,4);
3358 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
3364 volume+=thickness*gauss->
weight*Jdet;
3368 xDelete<IssmDouble>(xyz_list);
3369 xDelete<IssmDouble>(values);
3371 return rho_ice*volume;
3385 if (segment_id!=this->
id)
_error_(
"error message: segment with id " << segment_id <<
" does not belong to element with id:" << this->
id);
3399 IssmDouble length=sqrt(pow(x2-x1,2)+pow(y2-y1,2));
3422 mass_flux= rho_ice*length*(
3423 (1./3.*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*nx+
3424 (1./3.*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*ny
3434 return this->
MassFlux(segment[0],segment[1],segment[2],segment[3],reCast<int>(segment[4]));
3457 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
3470 Jelem+=0.5*(model-observation)*(model-observation)*Jdet*weight*gauss->
weight;
3496 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
3507 Jelem+=Jdet*weight*gauss->
weight;
3527 return new GaussTria(area_coordinates,order);
3532 return new GaussTria(point1,fraction1,fraction2,mainlyfloating,order);
3537 return new GaussTria(point1,fraction1,fraction2,order);
3544 return new GaussTria(area_coordinates,order_vert);
3551 return new GaussTria(indices[0],indices[1],order);
3558 return new GaussTria(indices[0],indices[1],order);
3651 if(found)*pvalue=value;
3661 vector[0]=xyz_list[1*3+0] - xyz_list[0*3+0];
3662 vector[1]=xyz_list[1*3+1] - xyz_list[0*3+1];
3664 norm=sqrt(vector[0]*vector[0] + vector[1]*vector[1]);
3666 bed_normal[0]= + vector[1]/norm;
3667 bed_normal[1]= - vector[0]/norm;
3677 vector[0]=xyz_list[1*3+0] - xyz_list[0*3+0];
3678 vector[1]=xyz_list[1*3+1] - xyz_list[0*3+1];
3680 norm=sqrt(vector[0]*vector[0] + vector[1]*vector[1]);
3682 normal[0]= + vector[1]/norm;
3683 normal[1]= - vector[0]/norm;
3694 vector[0]=xyz_list[1*3+0] - xyz_list[0*3+0];
3695 vector[1]=xyz_list[1*3+1] - xyz_list[0*3+1];
3697 norm=sqrt(vector[0]*vector[0] + vector[1]*vector[1]);
3699 top_normal[0]= + vector[1]/norm;
3700 top_normal[1]= - vector[0]/norm;
3727 density=rho_ice/rho_water;
3736 bed_hydro=-density*h[i];
3754 int indices[2]={6,7};
3761 int* indices=xNew<int>(size);
3762 for(
int i=0;i<size;i++) indices[i] = offset+i;
3764 xDelete<int>(indices);
3770 int indices[2]={6,7};
3777 int* indices=xNew<int>(size);
3778 for(
int i=0;i<size;i++) indices[i] = offset+i;
3780 xDelete<int>(indices);
3802 vertexonbase = xNew<IssmDouble>(numnodes);
3809 if(vertexonbase[i]==1){
3812 groundedicelevelset_input->
GetInputValue(&groundedice,gauss);
3816 xz_plane[0]=cos(theta); xz_plane[3]=0.;
3817 xz_plane[1]=sin(theta); xz_plane[4]=0.;
3818 xz_plane[2]=0.; xz_plane[5]=1.;
3832 xDelete<IssmDouble>(vertexonbase);
3888 if(basin_icefront_area[reCast<int>(basinid[iv])-1]==0.) meltrates[iv]=0.;
3891 qsg_basin=qsg/basin_icefront_area[reCast<int>(basinid[iv])-1];
3894 meltrates[iv]=((A*
max(-bed,0.)*pow(
max(qsg_basin,0.),
alpha)+B)*pow(
max(TF,0.),beta))/86400;
3897 if(xIsNan<IssmDouble>(meltrates[iv]))
_error_(
"NaN found in vector");
3898 if(xIsInf<IssmDouble>(meltrates[iv]))
_error_(
"Inf found in vector");
3905 xDelete<IssmDouble>(basin_icefront_area);
3920 control_init=control_enum;
3939 for(
int n=0;n<N;n++){
3942 values[i]=vector[idlist[i]];
3953 else _error_(
"Type not supported");
3988 values[i]=vector[idlist[i]];
3999 int analysis_counter;
4046 _error_(
"not implemented yet");
4052 int analysis_counter;
4065 this->
SpawnSegHook(xDynamicCast<ElementHook*>(seg),index1,index2);
4096 _error_(
"not implemented yet");
4129 this->
StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
4130 strainxx=epsilon[0];
4131 strainyy=epsilon[1];
4132 strainxy=epsilon[2];
4135 strainparallel[iv]=(vx*vx*(strainxx)+vy*vy*(strainyy)+2*vy*vx*strainxy)/(vel+1.e-14);
4143 xDelete<IssmDouble>(xyz_list);
4175 this->
StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
4176 strainxx=epsilon[0];
4177 strainyy=epsilon[1];
4178 strainxy=epsilon[2];
4181 strainperpendicular[iv]=(vx*vx*(strainyy)+vy*vy*(strainxx)-2*vy*vx*strainxy)/(vel+1.e-14);
4189 xDelete<IssmDouble>(xyz_list);
4204 for(
int i=0;i<3;i++){
4205 v13[i]=xyz_list[0][i]-xyz_list[2][i];
4206 v23[i]=xyz_list[1][i]-xyz_list[2][i];
4209 normal[0]=v13[1]*v23[2]-v13[2]*v23[1];
4210 normal[1]=v13[2]*v23[0]-v13[0]*v23[2];
4211 normal[2]=v13[0]*v23[1]-v13[1]*v23[0];
4213 S = 0.5 * sqrt(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]);
4247 if(xyz_list[i][0]<minx) minx=xyz_list[i][0];
4248 if(xyz_list[i][0]>maxx) maxx=xyz_list[i][0];
4249 if(xyz_list[i][1]<miny) miny=xyz_list[i][1];
4250 if(xyz_list[i][1]>maxy) maxy=xyz_list[i][1];
4269 int domaintype,index1,index2;
4284 if(gl[0]==0.) gl[0]=gl[0]+epsilon;
4285 if(gl[1]==0.) gl[1]=gl[1]+epsilon;
4286 if(gl[2]==0.) gl[2]=gl[2]+epsilon;
4297 s1=gl[2]/(gl[2]-gl[1]);
4298 s2=gl[2]/(gl[2]-gl[0]);
4302 xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
4303 xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
4304 xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
4305 xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
4306 xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
4307 xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
4309 else if(gl[1]*gl[2]>0){
4312 s1=gl[0]/(gl[0]-gl[1]);
4313 s2=gl[0]/(gl[0]-gl[2]);
4318 xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
4319 xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
4320 xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
4321 xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
4322 xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
4323 xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
4325 else if(gl[0]*gl[2]>0){
4328 s1=gl[1]/(gl[1]-gl[0]);
4329 s2=gl[1]/(gl[1]-gl[2]);
4334 xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
4335 xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
4336 xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
4337 xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
4338 xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
4339 xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
4355 normal[0] = -normal[0];
4356 normal[1] = -normal[1];
4360 IssmDouble calvingratex,calvingratey,thickness,Jdet;
4363 Input2* calvingratex_input=NULL;
4364 Input2* calvingratey_input=NULL;
4376 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
4384 flux += rho_ice*Jdet*gauss->
weight*thickness*(calvingratex*normal[0] + calvingratey*normal[1]);
4398 int domaintype,index1,index2;
4413 if(gl[0]==0.) gl[0]=gl[0]+epsilon;
4414 if(gl[1]==0.) gl[1]=gl[1]+epsilon;
4415 if(gl[2]==0.) gl[2]=gl[2]+epsilon;
4426 s1=gl[2]/(gl[2]-gl[1]);
4427 s2=gl[2]/(gl[2]-gl[0]);
4431 xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
4432 xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
4433 xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
4434 xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
4435 xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
4436 xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
4438 else if(gl[1]*gl[2]>0){
4441 s1=gl[0]/(gl[0]-gl[1]);
4442 s2=gl[0]/(gl[0]-gl[2]);
4447 xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
4448 xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
4449 xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
4450 xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
4451 xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
4452 xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);
4454 else if(gl[0]*gl[2]>0){
4457 s1=gl[1]/(gl[1]-gl[0]);
4458 s2=gl[1]/(gl[1]-gl[2]);
4463 xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
4464 xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
4465 xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
4466 xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
4467 xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
4468 xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
4484 normal[0] = -normal[0];
4485 normal[1] = -normal[1];
4489 IssmDouble calvingratex,calvingratey,vx,vy,vel,meltingrate,meltingratex,meltingratey,thickness,Jdet;
4492 Input2* calvingratex_input=NULL;
4493 Input2* calvingratey_input=NULL;
4496 Input2* meltingrate_input=NULL;
4511 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
4521 meltingratex=meltingrate*vx/(sqrt(vel)+1.e-14);
4522 meltingratey=meltingrate*vy/(sqrt(vel)+1.e-14);
4525 flux += rho_ice*Jdet*gauss->
weight*thickness*((calvingratex+meltingratex)*normal[0] + (calvingratey+meltingratey)*normal[1]);
4535 bool mainlyfloating;
4537 IssmDouble rho_ice,fraction1,fraction2,floatingmelt,Jdet,scalefactor;
4540 Gauss* gauss = NULL;
4548 Input2* scalefactor_input = NULL;
4556 gauss = this->
NewGauss(point1,fraction1,fraction2,1-mainlyfloating,3);
4557 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
4566 fbmb+=floatingmelt*Jdet*gauss->
weight*scalefactor;
4569 Total_Fbmb=rho_ice*fbmb;
4580 bool mainlyfloating;
4582 IssmDouble rho_ice,fraction1,fraction2,groundedmelt,Jdet,scalefactor;
4585 Gauss* gauss = NULL;
4593 Input2* scalefactor_input = NULL;
4601 gauss = this->
NewGauss(point1,fraction1,fraction2,mainlyfloating,2);
4602 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
4611 gbmb+=groundedmelt*Jdet*gauss->
weight*scalefactor;
4614 Total_Gbmb=rho_ice*gbmb;
4638 base = 1./2. * fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
4650 Total_Smb=rho_ice*base*smb*scalefactor;
4660 int* tria_node_ids = NULL;
4670 switch(finiteelement_type){
4673 tria_node_ids = xNew<int>(numnodes);
4674 tria_node_ids[0]= index + 1;
4678 tria_node_ids = xNew<int>(numnodes);
4679 tria_node_ids[0]=iomodel->
elements[3*index+0];
4680 tria_node_ids[1]=iomodel->
elements[3*index+1];
4681 tria_node_ids[2]=iomodel->
elements[3*index+2];
4685 tria_node_ids = xNew<int>(numnodes);
4686 tria_node_ids[0]=3*index+1;
4687 tria_node_ids[1]=3*index+2;
4688 tria_node_ids[2]=3*index+3;
4692 tria_node_ids = xNew<int>(numnodes);
4693 tria_node_ids[0]=iomodel->
elements[3*index+0];
4694 tria_node_ids[1]=iomodel->
elements[3*index+1];
4695 tria_node_ids[2]=iomodel->
elements[3*index+2];
4700 tria_node_ids = xNew<int>(numnodes);
4701 tria_node_ids[0]=iomodel->
elements[3*index+0];
4702 tria_node_ids[1]=iomodel->
elements[3*index+1];
4703 tria_node_ids[2]=iomodel->
elements[3*index+2];
4710 tria_node_ids = xNew<int>(numnodes);
4711 tria_node_ids[0]=iomodel->
elements[3*index+0];
4712 tria_node_ids[1]=iomodel->
elements[3*index+1];
4713 tria_node_ids[2]=iomodel->
elements[3*index+2];
4721 tria_node_ids = xNew<int>(numnodes);
4722 tria_node_ids[0]=iomodel->
elements[3*index+0];
4723 tria_node_ids[1]=iomodel->
elements[3*index+1];
4724 tria_node_ids[2]=iomodel->
elements[3*index+2];
4732 tria_node_ids = xNew<int>(numnodes);
4733 tria_node_ids[0]=iomodel->
elements[3*index+0];
4734 tria_node_ids[1]=iomodel->
elements[3*index+1];
4735 tria_node_ids[2]=iomodel->
elements[3*index+2];
4745 tria_node_ids = xNew<int>(numnodes);
4746 tria_node_ids[0]=iomodel->
elements[3*index+0];
4747 tria_node_ids[1]=iomodel->
elements[3*index+1];
4748 tria_node_ids[2]=iomodel->
elements[3*index+2];
4759 tria_node_ids = xNew<int>(numnodes);
4760 tria_node_ids[0]=iomodel->
elements[3*index+0];
4761 tria_node_ids[1]=iomodel->
elements[3*index+1];
4762 tria_node_ids[2]=iomodel->
elements[3*index+2];
4769 tria_node_ids = xNew<int>(numnodes);
4770 tria_node_ids[0]=iomodel->
elements[3*index+0];
4771 tria_node_ids[1]=iomodel->
elements[3*index+1];
4772 tria_node_ids[2]=iomodel->
elements[3*index+2];
4784 tria_node_ids = xNew<int>(numnodes);
4785 tria_node_ids[0]=iomodel->
elements[3*index+0];
4786 tria_node_ids[1]=iomodel->
elements[3*index+1];
4787 tria_node_ids[2]=iomodel->
elements[3*index+2];
4799 xDelete<int>(tria_node_ids);
4856 if (reCast<bool>(vertices_potentially_ungrounding[
vertices[i]->Pid()])){
4860 if(nodes_on_iceshelf[
vertices[i]->Pid()]>=0.){
4897 if(lsf[i]>maxvalue) maxvalue = lsf[i];
4898 if(lsf[i]<minvalue) minvalue = lsf[i];
4900 if(minvalue>fieldvalue)
return;
4901 if(maxvalue<fieldvalue)
return;
4906 int *indices = NULL;
4915 for(
int i=0;i<numiceverts;i++){
4917 x[counter] = xyz_list[indices[i]][0]+s[counter]*(xyz_list[indices[n]][0]-xyz_list[indices[i]][0]);
4918 y[counter] = xyz_list[indices[i]][1]+s[counter]*(xyz_list[indices[n]][1]-xyz_list[indices[i]][1]);
4926 if(lsf[indices[i]]==0.){
4927 x[counter]=xyz_list[indices[i]][0];
4928 y[counter]=xyz_list[indices[i]][1];
4931 if(counter==2)
break;
4941 _error_(
"not sure what's going on here...");
4949 xDelete<int>(indices);
4959 IssmDouble lithosphere_thickness,mantle_viscosity;
4975 int cross_section_shape;
4997 if (!mantle_viscosity_input)
_error_(
"mantle viscosity input needed to compute gia deflection!");
5002 if (!lithosphere_thickness_input)
_error_(
"lithosphere thickness input needed to compute gia deflection!");
5013 IssmDouble x0=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0;
5014 IssmDouble y0=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0;
5019 arguments.
times=times;
5030 arguments.
iedge=cross_section_shape;
5033 for(
int i=0;i<gsize;i++){
5037 IssmDouble ri=sqrt(pow(xi-x0,2)+pow(yi-y0,2));
5051 xDelete<IssmDouble>(hes);
5052 xDelete<IssmDouble>(times);
5082 bool store_green_functions=
false;
5086 if (!deltathickness_input)
_error_(
"delta thickness input needed to compute elastic adjustment!");
5108 x_element=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0;
5109 y_element=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0;
5118 U_elastic=xNewZeroInit<IssmDouble>(gsize);
5119 N_elastic=xNewZeroInit<IssmDouble>(gsize);
5120 E_elastic=xNewZeroInit<IssmDouble>(gsize);
5121 X_elastic=xNewZeroInit<IssmDouble>(gsize);
5122 Y_elastic=xNewZeroInit<IssmDouble>(gsize);
5124 int* indices=xNew<int>(gsize);
5125 IssmDouble* U_values=xNewZeroInit<IssmDouble>(gsize);
5126 IssmDouble* N_values=xNewZeroInit<IssmDouble>(gsize);
5127 IssmDouble* E_values=xNewZeroInit<IssmDouble>(gsize);
5128 IssmDouble* X_values=xNewZeroInit<IssmDouble>(gsize);
5129 IssmDouble* Y_values=xNewZeroInit<IssmDouble>(gsize);
5133 for(
int i=0;i<gsize;i++){
5141 dx = x_element - xx[i]; dy = y_element - yy[i];
5142 dist = sqrt(pow(dx,2)+pow(dy,2));
5143 alpha = dist*360.0/(2*
PI*earth_radius) *
PI/180.0;
5146 ang =
PI/2 - atan2(dy,dx);
5151 int index=reCast<int,IssmDouble>(
alpha/
PI*(M-1));
5152 U_elastic[i] += U_elastic_precomputed[index];
5153 Y_elastic[i] += H_elastic_precomputed[index]*Y_azim;
5154 X_elastic[i] += H_elastic_precomputed[index]*X_azim;
5157 U_values[i]+=3*rho_ice/rho_earth*area/(4*
PI*pow(earth_radius,2))*I*U_elastic[i];
5158 Y_values[i]+=3*rho_ice/rho_earth*area/(4*
PI*pow(earth_radius,2))*I*Y_elastic[i];
5159 X_values[i]+=3*rho_ice/rho_earth*area/(4*
PI*pow(earth_radius,2))*I*X_elastic[i];
5163 ang2 =
PI/2 - atan2(yy[i],xx[i]);
5165 else if (hemi == 1) {
5166 ang2 =
PI/2 - atan2(-yy[i],-xx[i]);
5169 N_azim = Y_azim*cos(ang2) + X_azim*sin(ang2);
5170 E_azim = X_azim*cos(ang2) - Y_azim*sin(ang2);
5171 N_elastic[i] += H_elastic_precomputed[index]*N_azim;
5172 E_elastic[i] += H_elastic_precomputed[index]*E_azim;
5173 N_values[i]+=3*rho_ice/rho_earth*area/(4*
PI*pow(earth_radius,2))*I*N_elastic[i];
5174 E_values[i]+=3*rho_ice/rho_earth*area/(4*
PI*pow(earth_radius,2))*I*E_elastic[i];
5185 xDelete<int>(indices);
5186 xDelete<IssmDouble>(U_values); xDelete<IssmDouble>(N_values); xDelete<IssmDouble>(E_values);
5187 xDelete<IssmDouble>(U_elastic); xDelete<IssmDouble>(N_elastic); xDelete<IssmDouble>(E_elastic);
5188 xDelete<IssmDouble>(X_values); xDelete<IssmDouble>(Y_values);
5189 xDelete<IssmDouble>(X_elastic); xDelete<IssmDouble>(Y_elastic);
5198 bool spherical=
true;
5220 bool store_green_functions=
false;
5224 if (!deltathickness_input)
_error_(
"delta thickness input needed to compute elastic adjustment!");
5247 minlong=400; maxlong=-20;
5249 llr_list[i][0]=(90-llr_list[i][0]);
5250 if(llr_list[i][1]<0)llr_list[i][1]=180+(180+llr_list[i][1]);
5251 if(llr_list[i][1]>maxlong)maxlong=llr_list[i][1];
5252 if(llr_list[i][1]<minlong)minlong=llr_list[i][1];
5254 if(minlong==0 && maxlong>180){
5255 if (llr_list[0][1]==0)llr_list[0][1]=360;
5256 if (llr_list[1][1]==0)llr_list[1][1]=360;
5257 if (llr_list[2][1]==0)llr_list[2][1]=360;
5262 if(llr_list[0][0]==0)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
5263 if(llr_list[1][0]==0)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
5264 if(llr_list[2][0]==0)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;
5268 if(llr_list[0][0]==180)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
5269 if(llr_list[1][0]==180)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
5270 if(llr_list[2][0]==180)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;
5272 late=(llr_list[0][0]+llr_list[1][0]+llr_list[2][0])/3.0;
5273 longe=(llr_list[0][1]+llr_list[1][1]+llr_list[2][1])/3.0;
5276 if(longe>180)longe=longe-360;
5285 x_element=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0;
5286 y_element=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0;
5287 z_element=(xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2])/3.0;
5296 U_elastic=xNewZeroInit<IssmDouble>(gsize);
5297 N_elastic=xNewZeroInit<IssmDouble>(gsize);
5298 E_elastic=xNewZeroInit<IssmDouble>(gsize);
5300 int* indices=xNew<int>(gsize);
5301 IssmDouble* U_values=xNewZeroInit<IssmDouble>(gsize);
5302 IssmDouble* N_values=xNewZeroInit<IssmDouble>(gsize);
5303 IssmDouble* E_values=xNewZeroInit<IssmDouble>(gsize);
5309 for(
int i=0;i<gsize;i++){
5314 lati=latitude[i]/180*
PI; longi=longitude[i]/180*
PI;
5316 delPhi=fabs(lati-late); delLambda=fabs(longi-longe);
5317 alpha=2.*asin(sqrt(pow(sin(delPhi/2),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
5320 x = xx[i]; y = yy[i]; z = zz[i];
5321 if(latitude[i]==90){
5324 if(latitude[i]==-90){
5327 dx = x_element-x; dy = y_element-y; dz = z_element-z;
5328 N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
5329 E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
5332 int index=reCast<int,IssmDouble>(
alpha/
PI*(M-1));
5333 U_elastic[i] += U_elastic_precomputed[index];
5334 N_elastic[i] += H_elastic_precomputed[index]*N_azim;
5335 E_elastic[i] += H_elastic_precomputed[index]*E_azim;
5338 U_values[i]+=3*rho_ice/rho_earth*area/planetarea*I*U_elastic[i];
5339 N_values[i]+=3*rho_ice/rho_earth*area/planetarea*I*N_elastic[i];
5340 E_values[i]+=3*rho_ice/rho_earth*area/planetarea*I*E_elastic[i];
5347 xDelete<int>(indices);
5348 xDelete<IssmDouble>(U_values); xDelete<IssmDouble>(N_values); xDelete<IssmDouble>(E_values);
5349 xDelete<IssmDouble>(U_elastic); xDelete<IssmDouble>(N_elastic); xDelete<IssmDouble>(E_elastic);
5355 #ifdef _HAVE_SEALEVELRISE_
5392 bool spherical=
true;
5401 llr_list[i][0]=(90-llr_list[i][0]);
5402 if(llr_list[i][1]<0)llr_list[i][1]=180+(180+llr_list[i][1]);
5403 if(llr_list[i][1]>maxlong)maxlong=llr_list[i][1];
5404 if(llr_list[i][1]<minlong)minlong=llr_list[i][1];
5406 if(minlong==0 && maxlong>180){
5407 if (llr_list[0][1]==0)llr_list[0][1]=360;
5408 if (llr_list[1][1]==0)llr_list[1][1]=360;
5409 if (llr_list[2][1]==0)llr_list[2][1]=360;
5413 if(llr_list[0][0]==0)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
5414 if(llr_list[1][0]==0)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
5415 if(llr_list[2][0]==0)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;
5418 if(llr_list[0][0]==180)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
5419 if(llr_list[1][0]==180)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
5420 if(llr_list[2][0]==180)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;
5422 late=(llr_list[0][0]+llr_list[1][0]+llr_list[2][0])/3.0;
5423 longe=(llr_list[0][1]+llr_list[1][1]+llr_list[2][1])/3.0;
5426 if(longe>180)longe=(longe-180)-180;
5431 re=(llr_list[0][2]+llr_list[1][2]+llr_list[2][2])/3.0;
5447 dI_list[0] = -4*
PI*(rho_water*S*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/planetarea;
5448 dI_list[1] = -4*
PI*(rho_water*S*area)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/planetarea;
5449 dI_list[2] = +4*
PI*(rho_water*S*area)*pow(re,4)*(1-pow(sin(late),2))/planetarea;
5459 if (!deltathickness_input)
_error_(
"delta thickness input needed to compute sea level rise!");
5462 dI_list[0] = -4*
PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/planetarea;
5463 dI_list[1] = -4*
PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/planetarea;
5464 dI_list[2] = +4*
PI*(rho_ice*I*area)*pow(re,4)*(1-pow(sin(late),2))/planetarea;
5488 bool spherical=
true;
5497 IssmDouble x_element,y_element,z_element,x,y,z,dx,dy,dz,N_azim,E_azim;
5499 #ifdef _HAVE_RESTRICT_
5504 IssmDouble* __restrict__ G_elastic_precomputed=NULL;
5505 IssmDouble* __restrict__ G_rigid_precomputed=NULL;
5506 IssmDouble* __restrict__ U_elastic_precomputed=NULL;
5507 IssmDouble* __restrict__ H_elastic_precomputed=NULL;
5526 bool computerigid =
true;
5527 bool computeelastic =
true;
5553 indices=xNew<IssmDouble>(gsize);
5554 G=xNewZeroInit<IssmDouble>(gsize);
5555 GU=xNewZeroInit<IssmDouble>(gsize);
5557 GN=xNewZeroInit<IssmDouble>(gsize);
5558 GE=xNewZeroInit<IssmDouble>(gsize);
5571 x_element=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0;
5572 y_element=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0;
5573 z_element=(xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2])/3.0;
5576 late= asin(z_element/planetradius);
5577 longe = atan2(y_element,x_element);
5579 constant=3/rho_earth*area/planetarea;
5581 for(
int i=0;i<gsize;i++){
5586 lati=latitude[i]/180*
PI; longi=longitude[i]/180*
PI;
5587 delPhi=fabs(lati-late); delLambda=fabs(longi-longe);
if (delLambda>
PI)delLambda=2*
PI-delLambda;
5588 alpha=2.*asin(sqrt(pow(sin(delPhi/2),2)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
5589 indices[i]=
alpha/
PI*reCast<IssmDouble,int>(M-1);
5590 index=reCast<int,IssmDouble>(indices[i]);
5594 G[i]+=G_rigid_precomputed[index];
5597 G[i]+=G_elastic_precomputed[index];
5602 GU[i] = constant * U_elastic_precomputed[index];
5605 x = xx[i]; y = yy[i]; z = zz[i];
5606 if(latitude[i]==90){
5609 if(latitude[i]==-90){
5612 dx = x_element-x; dy = y_element-y; dz = z_element-z;
5613 N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
5614 E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
5616 GN[i] = constant*H_elastic_precomputed[index]*N_azim;
5617 GE[i] = constant*H_elastic_precomputed[index]*E_azim;
5632 #ifdef _HAVE_RESTRICT_
5655 int bp_compute_fingerprints= 0;
5659 if(bp_compute_fingerprints)this->SealevelriseBottomPressure(Sgi,masks);
5662 this->SealevelriseEustaticIce(Sgi,peustatic,masks, oceanarea);
5665 this->SealevelriseEustaticHydro(Sgi,peustatic,masks, oceanarea);
5677 bool notfullygrounded=
false;
5678 bool scaleoceanarea=
false;
5734 if(notfullygrounded){
5744 if (!deltathickness_input)
_error_(
"delta thickness input needed to compute sea level rise!");
5750 bool mainlyfloating =
true;
5756 Gauss* gauss = this->
NewGauss(point1,fraction1,fraction2,mainlyfloating,2);
5761 for(
int ig=gauss->
begin();ig<gauss->end();ig++){
5766 total_weight+=gauss->
weight;
5775 if(scaleoceanarea) oceanarea=3.619e+14;
5776 eustatic += rho_ice*area*phi*I/(oceanarea*rho_water);
5782 for(
int i=0;i<gsize;i++) Sgi[i]+=G[i]*I;
5785 _assert_(!xIsNan<IssmDouble>(eustatic));
5786 *peustatic=eustatic;
5797 bool notfullygrounded=
false;
5798 bool scaleoceanarea=
false;
5814 if(masks->
isiceonly[this->lid]){ *peustatic=0;
return; }
5817 if (masks->
isfullyfloating[this->lid]){ constant=0; *peustatic=0;
return; }
5836 if (!deltathickness_input)
_error_(
"SurfaceloadWaterHeightChangeEnum input needed to compute sea level rise!");
5841 if(scaleoceanarea) oceanarea=3.619e+14;
5842 eustatic += rho_freshwater*area*phi*W/(oceanarea*rho_water);
5845 W=W*rho_freshwater*phi;
5848 for(
int i=0;i<gsize;i++) Sgi[i]+=G[i]*W;
5851 _assert_(!xIsNan<IssmDouble>(eustatic));
5852 *peustatic=eustatic;
5889 if (!bottompressure_change_input)
_error_(
"bottom pressure input needed to compute sea level rise fingerprint!");
5896 for(
int i=0;i<gsize;i++) Sgi[i]+=G[i]*BP;
5909 int bp_compute_fingerprints= 0;
5940 for(
int i=0;i<gsize;i++) Sgo[i]+=G[i]*S;
5952 int bp_compute_fingerprints= 0;
5960 bool computeelastic=
true;
5971 if(!computeelastic)
return;
5992 for(
int i=0;i<gsize;i++){
6004 if (!deltathickness_input)
_error_(
"delta thickness input needed to compute sea level rise!");
6010 for(
int i=0;i<gsize;i++){
6024 #ifdef _HAVE_DAKOTA_
6025 void Tria::InputUpdateFromMatrixDakota(
IssmDouble* matrix,
int nrows,
int ncols,
int name,
int type){
6041 for(
int t=0;t<ncols;t++){
6042 for(
int i=0;i<3;i++){
6044 values[i]=matrix[ncols*row+t];
6056 for(
int t=0;t<ncols;t++){
6068 void Tria::InputUpdateFromVectorDakota(
IssmDouble* vector,
int name,
int type){
6108 di = rho_ice/rho_water;
6111 for (j=0; j<3; j++) {
6113 if (hydrostatic_ratio[j] >= 0.)
6114 thickness[j]=values[j]-(values[j]/thickness_init[j]-1.)*hydrostatic_ratio[j]*surface[j]/(1.-di);
6117 thickness[j]=thickness_init[j];
6120 if (thickness[j] < 0.) thickness[j]=1./(1.-di);
6121 bed[j]=surface[j]-thickness[j];
6126 for (j=0; j<3; j++) {
6128 if (hydrostatic_ratio[j] >= 0.)
6129 thickness[j]=values[j];
6132 thickness[j]=thickness_init[j];
6136 for(j=0;j<3;j++)bed[j]=surface[j]-thickness[j];
6156 value=vector[this->
Sid()];