Skip to content

Commit

Permalink
CHG: adding factors to analysis
Browse files Browse the repository at this point in the history
  • Loading branch information
@hseroussi
hseroussi committed Jun 24, 2024
1 parent 7885d2b commit 44d1c2c
Show file tree
Hide file tree
Showing 4 changed files with 30 additions and 20 deletions.
6 changes: 4 additions & 2 deletions src/c/analyses/SmoothAnalysis.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -85,9 +85,10 @@ ElementMatrix* SmoothAnalysis::CreateKMatrix(Element* element){/*{{{*/
element->NodalFunctions(basis,gauss);
element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);

IssmDouble factor = gauss->weight*Jdet;
for(int i=0;i<numnodes;i++){
for(int j=0;j<numnodes;j++){
Ke->values[i*numnodes+j] += gauss->weight*Jdet*(
Ke->values[i*numnodes+j] += factor*(
basis[i]*basis[j]
+(l*thickness)*(l*thickness)*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j])
);
Expand Down Expand Up @@ -194,7 +195,8 @@ ElementVector* SmoothAnalysis::CreatePVector(Element* element){/*{{{*/
input->GetInputValue(&value,gauss);
}

for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*value*basis[i];
IssmDouble factor = Jdet*gauss->weight*value;
for(int i=0;i<numnodes;i++) pe->values[i]+=factor*basis[i];
}

/*Clean up and return*/
Expand Down
9 changes: 5 additions & 4 deletions src/c/analyses/StressbalanceSIAAnalysis.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -443,13 +443,14 @@ ElementVector* StressbalanceSIAAnalysis::CreatePVector3D(Element* element){/*{{{
z = element->GetZcoord(xyz_list,gauss);
element->JacobianDeterminantLine(&Jdet,&xyz_list_line[0][0],gauss);

IssmDouble factor = constant_part*pow((surface-z)/B,n)*Jdet*gauss->weight/connectivity[1];
if(element->IsOnSurface()){
pe->values[2*nodeup+0]+=constant_part*pow((surface-z)/B,n)*slope[0]*Jdet*gauss->weight/connectivity[1];
pe->values[2*nodeup+1]+=constant_part*pow((surface-z)/B,n)*slope[1]*Jdet*gauss->weight/connectivity[1];
pe->values[2*nodeup+0]+=factor*slope[0];
pe->values[2*nodeup+1]+=factor*slope[1];
}
else{/*connectivity is too large, should take only half on it*/
pe->values[2*nodeup+0]+=constant_part*pow((surface-z)/B,n)*slope[0]*Jdet*gauss->weight*2./connectivity[1];
pe->values[2*nodeup+1]+=constant_part*pow((surface-z)/B,n)*slope[1]*Jdet*gauss->weight*2./connectivity[1];
pe->values[2*nodeup+0]+=factor*slope[0]*2.;
pe->values[2*nodeup+1]+=factor*slope[1]*2.;
}
}

Expand Down
12 changes: 8 additions & 4 deletions src/c/analyses/StressbalanceVerticalAnalysis.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -314,9 +314,10 @@ ElementMatrix* StressbalanceVerticalAnalysis::CreateKMatrixVolume(Element* eleme
element->NodalFunctions(basis,gauss);
element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);

IssmDouble factor = gauss->weight*Jdet;
for(int i=0;i<numnodes;i++){
for(int j=0;j<numnodes;j++){
Ke->values[i*numnodes+j] += gauss->weight*Jdet*(
Ke->values[i*numnodes+j] += factor*(
basis[j]*dbasis[2*numnodes+i]
);
}
Expand Down Expand Up @@ -405,7 +406,8 @@ ElementVector* StressbalanceVerticalAnalysis::CreatePVectorBase(Element* element
element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
element->NodalFunctions(basis,gauss);

for(int i=0;i<numnodes;i++) pe->values[i]+=-Jdet*gauss->weight*(vx*dbdx+vy*dbdy-vz-basalmeltingvalue)*basis[i];
IssmDouble factor = -Jdet*gauss->weight*(vx*dbdx+vy*dbdy-vz-basalmeltingvalue);;
for(int i=0;i<numnodes;i++) pe->values[i]+=factor*basis[i];
}

/*Clean up and return*/
Expand Down Expand Up @@ -464,7 +466,8 @@ ElementVector* StressbalanceVerticalAnalysis::CreatePVectorSurface(Element* elem
element->JacobianDeterminantTop(&Jdet,xyz_list_surface,gauss);
element->NodalFunctions(basis,gauss);

for(int i=0;i<numnodes;i++) pe->values[i]+=-Jdet*gauss->weight*(vx*dsdx+vy*dsdy-vz+smb)*basis[i];
IssmDouble factor = -Jdet*gauss->weight*(vx*dsdx+vy*dsdy-vz+smb);
for(int i=0;i<numnodes;i++) pe->values[i]+=factor*basis[i];
}

/*Clean up and return*/
Expand Down Expand Up @@ -516,7 +519,8 @@ ElementVector* StressbalanceVerticalAnalysis::CreatePVectorVolume(Element* eleme
dudx=du[0];
dvdy=dv[1];

for(int i=0;i<numnodes;i++) pe->values[i] += (dudx+dvdy+dwdz)*Jdet*gauss->weight*basis[i];
IssmDouble factor = (dudx+dvdy+dwdz)*Jdet*gauss->weight;
for(int i=0;i<numnodes;i++) pe->values[i] += factor*basis[i];
}

/*Clean up and return*/
Expand Down
23 changes: 13 additions & 10 deletions src/c/analyses/ThermalAnalysis.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -294,7 +294,7 @@ ElementMatrix* ThermalAnalysis::CreateKMatrixVolume(Element* element){/*{{{*/
int stabilization;
IssmDouble Jdet,dt,u,v,w,um,vm,wm,vel;
IssmDouble h,hx,hy,hz,vx,vy,vz,D_scalar;
IssmDouble tau_parameter,diameter;
IssmDouble tau_parameter,diameter,factor;
IssmDouble tau_parameter_anisotropic[2],tau_parameter_hor,tau_parameter_ver;
IssmDouble* xyz_list = NULL;

Expand Down Expand Up @@ -336,11 +336,11 @@ ElementMatrix* ThermalAnalysis::CreateKMatrixVolume(Element* element){/*{{{*/
if(dt!=0.) D_scalar=D_scalar*dt;

/*Conduction: */
factor = D_scalar*kappa;
for(int i=0;i<numnodes;i++){
for(int j=0;j<numnodes;j++){
Ke->values[i*numnodes+j] += D_scalar*kappa*(
dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[2*numnodes+j]*dbasis[2*numnodes+i]
);
Ke->values[i*numnodes+j] += factor*(
dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[2*numnodes+j]*dbasis[2*numnodes+i]);
}
}

Expand Down Expand Up @@ -368,9 +368,10 @@ ElementMatrix* ThermalAnalysis::CreateKMatrixVolume(Element* element){/*{{{*/
element->ElementSizes(&hx,&hy,&hz);
vel=sqrt(vx*vx + vy*vy + vz*vz)+1.e-14;
h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) + pow(hz*vz/vel,2));
K[0][0]=h/(2.*vel)*fabs(vx*vx); K[0][1]=h/(2.*vel)*fabs(vx*vy); K[0][2]=h/(2.*vel)*fabs(vx*vz);
K[1][0]=h/(2.*vel)*fabs(vy*vx); K[1][1]=h/(2.*vel)*fabs(vy*vy); K[1][2]=h/(2.*vel)*fabs(vy*vz);
K[2][0]=h/(2.*vel)*fabs(vz*vx); K[2][1]=h/(2.*vel)*fabs(vz*vy); K[2][2]=h/(2.*vel)*fabs(vz*vz);
factor = h/(2.*vel);
K[0][0]=factor*fabs(vx*vx); K[0][1]=factor*fabs(vx*vy); K[0][2]=factor*fabs(vx*vz);
K[1][0]=factor*fabs(vy*vx); K[1][1]=factor*fabs(vy*vy); K[1][2]=factor*fabs(vy*vz);
K[2][0]=factor*fabs(vz*vx); K[2][1]=factor*fabs(vz*vy); K[2][2]=factor*fabs(vz*vz);
for(int i=0;i<3;i++) for(int j=0;j<3;j++) K[i][j] = D_scalar*K[i][j];

for(int i=0;i<numnodes;i++){
Expand All @@ -386,18 +387,20 @@ ElementMatrix* ThermalAnalysis::CreateKMatrixVolume(Element* element){/*{{{*/
else if(stabilization==2){
diameter=element->MinEdgeLength(xyz_list);
tau_parameter=element->StabilizationParameter(u-um,v-vm,w-wm,diameter,kappa);

factor = tau_parameter*D_scalar;
for(int i=0;i<numnodes;i++){
for(int j=0;j<numnodes;j++){
Ke->values[i*numnodes+j]+=tau_parameter*D_scalar*
Ke->values[i*numnodes+j]+=factor*
((u-um)*dbasis[0*numnodes+i]+(v-vm)*dbasis[1*numnodes+i]+(w-wm)*dbasis[2*numnodes+i])*
((u-um)*dbasis[0*numnodes+j]+(v-vm)*dbasis[1*numnodes+j]+(w-wm)*dbasis[2*numnodes+j]);
}
}
if(dt!=0.){
D_scalar=gauss->weight*Jdet;
D_scalar=gauss->weight*Jdet*tau_parameter;
for(int i=0;i<numnodes;i++){
for(int j=0;j<numnodes;j++){
Ke->values[i*numnodes+j]+=tau_parameter*D_scalar*basis[j]*((u-um)*dbasis[0*numnodes+i]+(v-vm)*dbasis[1*numnodes+i]+(w-wm)*dbasis[2*numnodes+i]);
Ke->values[i*numnodes+j]+=D_scalar*basis[j]*((u-um)*dbasis[0*numnodes+i]+(v-vm)*dbasis[1*numnodes+i]+(w-wm)*dbasis[2*numnodes+i]);
}
}
}
Expand Down

0 comments on commit 44d1c2c

Please sign in to comment.