Conversion from array to finite element space variable

General Discussion
1

Hi,

Variable initialization:

fespace Nh(Th, [P1, P1]);
// FEM variables.
Nh [u1, u2], [v1, v2];

So elasticity problem can be solved like this

problem elas([u1,u2],[v1,v2]) =
  int2d(Th)(  
	    lam*div(u1,u2)*div(v1,v2)	
	    +2.*mu*( epsilon(u1,u2)'*epsilon(v1,v2) )
	      )
  - int1d(Th, 2)(p * v1)
  + on(4, u1 = 0)
  + on(1, u2 = 0)
  ;

or equivalently using the variational definition like

  varf elas([u1, u2], [v1, v2]) =
    int2d(Th)(  
      lam * div(u1, u2) * div(v1, v2)	
      + 2. * mu * ( epsilon(u1, u2)' * epsilon(v1, v2)))
    + on(4, u1 = 0)
    + on(1, u2 = 0)
    ;

  varf l([u1, u2], [v1, v2]) = int1d(Th, 2)(loadFraction * p * v1)  
    + on(4, u1 = 0)
    + on(1, u2 = 0)
    ;
  matrix K = elas(Nh, Nh); // Stiffness matrix.
  real[int] b = l(0, Nh); // External forces and boundary conditions.
  real[int] U(K.n);
  U = K^-1 * b;

Problem

The first solution results in u1 and u2 being the displacements in x and y directions. While the second solution returns a large array with displacements in the form U = [ux1, uy1, .... , uxn, uyn].

Running macro over the displacements

macro e(u1, u2) [dx(u1), dy(u2), dy(u1) + dx(u2)] // EOM

works with the displacements u1 and u2, but not with the array U.

Question

How can I convert array U into types of u1 and u2? Is that possible?

2

If you want to get of u1 and u2 from the array U, do u1[]=U;. This copies U into the array degrees of freedom of [u1,u2] that belong to your fespace Nh.

3

Ehm, Maybe I am doing it wrong but this does not seem to work.

So U has some values. Which are OK. And then I run

macro e(u1, u2)
	[
		dx(u1),
		(dy(u1) + dx(u2)) / 2.,
		(dx(u2) + dy(u1)) / 2., 
		dy(u2)
	] // EOM
macro sigma(u1, u2)
	[
		(lam + 2. * G) * e(u1, u2)[0] + lam * e(u1, u2)[3],
		2. * G * e(u1, u2)[1],
		2. * G * e(u1, u2)[2],
		lam * e(u1, u2)[0] + (lam + 2. * G) * e(u1, u2)[3]
	] // EOM
  u1[] = U;

  // Stresses 
  fespace Stress(Th, P1);
  Stress sxx, syy, sxy;

  sxx = sigma(u1, u2)[0];
  syy = sigma(u1, u2)[3]; 
  sxy = sigma(u1, u2)[1];
  cout << "sxx" << endl;
  cout << sxx[].max << endl;
  cout << sxx[].min << endl;
  cout << "syy" << endl;
  cout << syy[].max << endl;
  cout << syy[].min << endl;
  cout << "sxy" << endl;
  cout << sxy[].max << endl;
  cout << sxy[].min << endl;

Which results in all values being 0, which is for sure not OK. Where the only step I am not sure about is the u1[] = U.

Also note

Also note that

u1[] = U

does not produce the same results as

  for(int i = 0; i < K.n / 2; i++){
      u1[](i) = U(2 * i);
      u2[](i) = U(2 * i + 1);
  }

So confusing…