How to apply the boundary conditons to the mass and stiffness matrices?

I'm using the Partial Differential equation toolbox for getting the mass and stiffness matrices of a cube and liked to apply boundary conditions on the faces. My code looks like that:

gm = multicuboid(2,2,2);
model = createpde;
model.Geometry = gm;


specifyCoefficients(model,'m',0,'d',1,'c',1,'a',0,'f',1);
applyBoundaryCondition(model,'dirichlet','Face',1:6,'u',0);
mesh = generateMesh(model,'GeometricOrder','linear', ...
    'Hmax', 2^-2);
FEM = asembleFEMatricess(model);   
M = FEM.M; K = FEM.K;

How can I apply the boundary conditions to the matrices M and K. I know, that FEM also contains the matrices

     FEM = 

  struct with fields:

    K: [919×919 double]
    A: [919×919 double]
    F: [919×1 double]
    Q: [919×919 double]
    G: [919×1 double]
    H: [452×919 double]
    R: [452×1 double]
    M: [919×919 double]

where G,H,R,M store some sort of information about the boundary conditions, but I'm not sure how to combine that with M and K. I'm also not sure if which matrices the 'nullspace' options returns.

PS: I need the matrices for solving the heat equation with a space-time-method, that is why a don't use the solve option from the model.

PPS: This should also work for quadratic meshs, thats why I can't just delete rows and columns.

[Sravan Kumar Putta]

Hello, even i am facing same issue. Did you resolve this problem. I would be great if you could help me

[Sravan Kumar Putta]

Once you got matrices from nullspace or stiff-spring, How did you solve it using space - time method ? If you have used ODE solvers then how did you frame the ode function?

Can any one explain me?

[Julian Henning]

I implemented a Crank-Nicolson scheme by hand for solving a wave equation. My code look like this:

NeN = length(edgeNodes); % number of edge nodes

% The starting condition (at t=0) is saved into U
U(NeN+1:end,1) = M\U0;

% The first step is calculated and saved into U
U(NeN+1:end,2) = U(NeN+1:end,1) + tau*(M\V0) + (tau*tau/2) * ...
    (M\((F_prev(NeN+1:end)-A*U(NeN+1:end,1))));

% The other time steps are calculated
for k=3:K
    % The calculated forces from the previous time step can be used again.
    % The previous force is now the previous previous force. The current
    % force is now the previous force. The current force has to be
    % calculated again. This process has the highest computational cost.
    
    F_prev_prev = F_prev;
    F_prev = F;
    
    if f_time_dependent
        F = zeros(N,1);
        for i=1:length(Elements)
            for j=1:3 % This inner loop yields better performance than an
                % equivalent vectorized expression
                F(Elements(j,i)) = F(Elements(j,i)) + 1/3 * F_area(i) *...
                    f(F_center(i,1),F_center(i,2),(k-1)*tau);
            end
        end
    end
        U(NeN+1:end, k) = (M+ tau^2 /4 * A)\...
            ((tau^2 / 4)* (F(NeN+1:end) + 2*F_prev(NeN+1:end) ...
            + F_prev_prev(NeN+1:end)) -(tau^2 / 4)*A*(2*U(NeN+1:end,k-1)...
            + U(NeN+1:end, k-2)) + M*(2*U(NeN+1:end,k-1)...
            - U(NeN+1:end, k-2)));
 
end

[Sravan Kumar Putta]

did you do it for heat equation?

[Sravan Kumar Putta]

I am in extreme need of help dear, can u pls look into this problem

https://in.mathworks.com/matlabcentral/answers/755649-how-to-solve-semi-discretized-pde-matrices-with-a-time-derivative-in-pde-tool-box-using-ode-solvers?s_tid=srchtitle

[Sravan Kumar Putta]

did you use nullspace or stiff-spring to impose the boundary condition?

Answer 1

[Ravi Kumar]

Use the 'nullspace' as second argument, you will get matrices with BC imposed by eleminating dirichlet DoFs.

Regards,

Ravi

[Julian Henning]

Thank you for your answer. This solution seems to work. I also tried 'stiff-spring' as an option, this also seems to work and for this option the size of the matrices still equals the number of nodes. If I use the 'nullspace' option the matrices have fewer elements (obviously the nodes with the Dirichlet conditions are missing). Is there an easy option to transfer the solution which also has fewer elements than the mesh back to the origianl size for easy plotting (meaning to add zeros at the 'correct' positions)?

Edit: I figured out that B * solutions seems to work. Is this correct?

[Ravi Kumar]

You have found the right way, B*u should expand u to full size.

[Sravan Kumar Putta]

Once you got matrices from nullspace or stiff-spring, How did you solve it using space - time method ? If you have used ODE solvers then how did you frame the ode function?

Can any one explain me?