Extracting data on boundary

General Discussion
1

Dear all,

I am solving a simple mechanical problem. I am able to solve the problem. I want to get the reactions at the fixed boundary. I tried but was not able to get the values correct. Can anyone tell me how to get reactions/stresses on the nodes?

// Solving simple mechanical problem

//Generating 2D mesh

int lower = 1;
int upper = 2;
int right = 3;
int left  = 4; 

border B(t=0,1){ x=0+t; y=0; label= lower; }
border S1(t=0,1){x= 1;y=0+t; label= right;}
border T(t=0,1){ x=1-t; y=1; label= upper; }
border S2(t=0,1){x= 0;y=1-t; label= left;}

int n=2;

plot(B(n) + S1(n) +T(n) +S2(n), wait=true);
mesh Omega = buildmesh(B(n) + S1(n) +T(n) +S2(n) );
plot(Omega,fill=1, cmm="Discritized Mesh",wait=true);
//---------------------------------------------------------------------------
// Finite element function space
fespace Displacement(Omega, P1); // linear shape functions 
Displacement ux, uy, delux, deluy;

fespace Stress(Omega, P0); // values to store
Stress sigmaxx, sigmayy, sigmaxy, Tx, Ty, Fx;      // To be checked

// Material Properties
real lambda = 121.15e3;
real mu = 80.77e3; 
real E= mu*(3*lambda+2*mu)/(lambda+mu);
real nuu = lambda/2/(lambda+mu);

// Constitutive realations
real sqrt2=sqrt(2.);
macro epsilon(ux,uy)  [dx(ux),dy(uy),(dy(ux)+dx(uy))/sqrt2] //EOM
macro div(ux,uy) (dx(ux)+dy(uy)) //EOM
macro C E/(1-nuu^2)*[[1,    nuu,                0],
                     [nuu,    1,                0],
                     [0,      0, 	   (1-nuu)/2]]  //EOM
macro Grad(s) [dx(s),dy(s)] //EOM

macro sigma(ux,uy)   lambda*div(ux,uy)*[1.0, 1.0, 0.0]+ 2*mu*epsilon(ux,uy) //EOM  

macro sigma1(ux,uy) [lambda*div(ux,uy)+2*mu*epsilon(ux,uy)[0],lambda*div(ux,uy)+2*mu*epsilon(ux,uy)[1], 
                     2*mu*epsilon(ux,uy)[0]] //EOM


// Applying loads

real ue=0.5;

problem mechanical([ux, uy],[delux, deluy]) =
        int2d(Omega)((C*epsilon(delux,deluy))'*epsilon(ux,uy)) 
      + on(left,ux=0,uy=0)
      + on(right,ux=ue)
;

//--------------------------------------------------------------------------------------
// solving problem

mechanical;

sigmaxx = sigma1(ux,uy)[0];
sigmayy = sigma1(ux,uy)[1];
sigmaxy = sigma1(ux,uy)[2];


plot(sigmaxx, cmm="Stresses", wait=true, fill=1);


Tx =0;
Ty =0;

Tx = int1d(Omega,left)(sigmaxx*N.x + sigmaxy*N.y);
Ty = sigmaxy*N.x + sigmayy*N.y;

plot(Tx, cmm="traction-x", wait=true, value=1, fill=1);
plot(Ty, cmm="traction-y", wait=true, value=1, fill=1);

Fx = int1d(Omega,left)(sigmaxx*N.x + sigmaxy*N.y);
//Fx = int1d(Omega,left)(sigmaxx*N.x + sigmaxy*N.y);

ofstream gg("data.dat");
gg << Tx[] <<endl;


real sum=0;

ofstream ff("values.dat");

for(int i=0;i<100;i++) {
	// x, y, Sxx, Syy, Sxy
		real yline = 0. + 1*i/100.;
		real xline = 1*i/100.;
        sum = sum + Fx(xline,yline);
	}

    ff << "Net reaction  ::" << sum <<endl;

plot(sigmaxx, cmm="Stresses", wait=true, fill=1);

cout<< "Solution Completed" << endl;
2

Lines 75 and 81 you assign a scalar value from int1d(...) to finite element space variables Tx and Fx. Is it intentional?
Furthermore you use P0 elements for stresses and forces, so they have constant value inside each element (only one degree of freedom).

1 Like