Definition of the c coefficient for assempde

Hello, I am working on an elasticity problem and I would like to solve the mechanical equilibrium equation with the pde toolbox (f=0, a=0). I work in 2D and have periodic boundary conditions. My definition of the inhomogeneous stiffness tensor (c in the pde) depends on a function defined in the whole domain (depending on x,y). The inhomogeneous stiffness tensor can be described through a 4x4 matrix (entries C_1111, C_1112, C1121...depending on the position). I have already tried to define the matrix as indicated in the manual, but also reducing my problem into a 3x3 matrix, my entries have dimensions of (nr. of triangles in the mesh (times) nr. of triangles in the mesh).. Can I use the pde toolbox to solve my problem? Thanks in advance!

Laura

Answer 1

[Bill Greene]

Hi,

Yes, you can solve your problem with PDE Toolbox.

I recommend writing a MATLAB function to define your c-coefficient as shown here for the scalar case:

http://www.mathworks.com/help/pde/ug/scalar-functional-form-and-calling-syntax.html

For your elasticity case, as you note, your c-coefficient can be thought of as a 4x4 matrix. There are several ways to define the terms in this matrix. If you want to define all the entries in the upper triangle, you would need to define 10 terms. As you note, the MATLAB function would need to return a matrix with dimensions 10 x number_of_triangles. Since your c-matrix is a function of x and y, each column in that matrix could have different values. The documentation page I list above shows how to calculate positions of the element centroids, solution values, and gradients that you might need to define your coefficients as a function of position.

Bill

[laureta]

Thanks for you answer, Bill! I have already tried to define it at each triangles' centroid in the mesh, defining the 10 terms [C11 0 C12], [0 C44;C44 0] and [C21 0 C22]. The thing is that each entry of the upper triangle has already dimensions of number_of_triangles x number_of_triangles because of the dependance of c on another function eta(x,y). I have a system of equations, so I need to define C11(x,y), C12(x,y)...right? . Can I call assempde(b,p,e,t,C11,C12,C21...) with the coefficients for the system of PDEs?

[Bill Greene]

Hi,

You can pass a 10 x number_of_triangles sized c-matrix to assempde like this assempde(b,p,e,t,c,...) Each of the terms in this matrix must be a scalar.

It sounds like you are running into difficulties computing the terms in your c-matrix from the results of your eta(x,y) function. If you post the code you use for this step, we may be able to suggest the appropriate code so you get a 10 x number_of_triangles matrix at the end of this step.

Bill

[laureta]

Hi Bill, My eta vector's components are defined as: eta1(1,1,1,1,1:n2,1:n1) = 0.5*(tanh((X-(L1/3))/0.5)-tanh((X-(2*L1/3))/0.5)); eta2(1,1,1,1,1:n2,1:n1) = 1 - eta1(1,1,1,1,1:n2,1:n1); where n1 and n2 are the xpts and ypts of the triangles of the mesh, L1 is the length of the domain, and [X,Y]=meshgrid(xpts,ypts). The dimensions are choosen to fit with the following loop: C1(i,j,k,l,:,:)=C1(i,j,k,l,:,:)+(eta1(1,1,1,1,:,:).^2).*A1(i,j,k,l) +(eta2(1,1,1,1,:,:).^2).*A2(i,j,k,l) ; A1 and A2 are scalars. As you can imagine, my stiffness tensor is C1. How can I reduce it to a 10 x nr_of_triangles matrix? Should I redefine anything? Thanks in advance!

Answer 2

[Bill Greene]

Hi,

Your coefficients are sufficiently complicated that I'm not sure I understand all the issues, but hopefully I can at least point you in a direction where you can resolve the details, yourself.

Looking at your code snippet, I see a couple of issues. Using 4 (or higher) dimensioned arrays probably simplifies your implementation from a mechanics point of view but presents an obstacle when translating to the form PDE Toolbox requires. Second, it is not useful to define the coefficients on a regular grid defined by meshgrid() since PDE Toolbox requires these coefficients at the centroids of each element in the model.

I've included a snippet of MATLAB code below that attempts to address these two issues. For testing, I just used matrices of ones for A1 and A2-- presumably you have the real values for these. You need to take a close look at this code to make sure I haven't made any coding errors and you may need to modify it for your specific case.

c can be defined this way and passed into assempde

A1 = ones(2,2,2,2); A2 = ones(2,2,2,2); % testing only!
c = @(p, t, u, time) calcCmat(p, t, A1, A2, L1);

The function calcCmat is defined by this code:

function c=calcCmat(p, t, A1, A2, L1)
% Triangle point indices
it1=t(1,:);
it2=t(2,:);
it3=t(3,:);

% Find centroids of triangles
X=(p(1,it1)+p(1,it2)+p(1,it3))/3;
eta1 = 0.5*(tanh((X-(L1/3))/0.5)-tanh((X-(2*L1/3))/0.5));
eta2 = 1 - eta1;
% map tensor components to PDE Toolbox c-matrix
ijkl = [1 1 1 1; 1 1 1 2; 1 1 2 2; 1 2 1 1; 1 2 2 1;
        2 2 1 1; 1 2 1 2; 1 2 2 2; 2 2 1 2; 2 2 2 2];
c = zeros(10, size(t,2));
for n=1:10
  i = ijkl(n,1); j = ijkl(n,2); k = ijkl(n,3); l = ijkl(n,4);
  c(n,:) = A1(i,j,k,l)*eta1.^2 + A2(i,j,k,l)*eta2.^2;
end
end

[laureta]

Thanks again for your advices! I have been trying to introduce your changes in my code. Unfortunately, I need to define all the 16 entries in the matrix because the stiffness matrix c is a sum for all entries. The original code was as follows, with Einstein summation:

         for i=1:2
            for j=1:2
                for k=1:2
                    for l=1:2
                        for p=1:2
                            for q=1:2
                                for r=1:2
                                    for s=1:2
                                            C1(i,j,k,l,:,:)=C1(i,j,k,l,:,:)+(eta1(1,1,1,1,:,:).^2).*(A1(i,p)*A1(j,q)*A1(k,r)*A1(l,s)*C(p,q,r,s,:,:)) ...
                                                +(eta2(1,1,1,1,:,:).^2).*(A2(i,p)*A2(j,q)*A2(k,r)*A2(l,s)*C(p,q,r,s,:,:)) ;
                                    end 
                                end 
                            end
                        end
                    end 
                end 
            end
        end  
where A1, A2 are rotation matrices. Sorry for not writing this sum from the beginning! I see it could have been easier... but thanks anyway! Do you recommend to keep trying to resolve it with this pde tool?

Answer 3

[Bill Greene]

>Do you recommend to keep trying to resolve it with this pde tool? Yes, definitely. I'm almost certain that the 10-entry version of the PDE Toolbox c-coefficient will support your case. It is simply a matter of mapping the entries of the 4'th order tensor, C1, to the ten entries in the matrix. I think that the mapping is defined by this vector ijkl = [1 1 1 1; 1 1 1 2; 1 1 2 2; 1 2 1 1; 1 2 2 1; 2 2 1 1; 1 2 1 2; 1 2 2 2; 2 2 1 2; 2 2 2 2]; where each row defines the indices for the 4th order tensor. But I could be wrong since I don't understand the details of your material. If your material is really strange and your c-tensor has no symmetry at all, you can use the 16-row version of the c-matrix which supports a completely general 4th order constitutive tensor in 2D.

Bill