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

Question

Julian Henning 2019-6-18

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/467638-how-to-apply-the-boundary-conditons-to-the-mass-and-stiffness-matrices

评论： Sravan Kumar Putta 2021-2-25

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.

6 个评论
显示 4更早的评论隐藏 4更早的评论

Julian Henning 2021-2-25

编辑：Julian Henning 2021-2-25

在 MATLAB Online 中打开

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 2021-2-25

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 2021-2-25

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

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Ravi Kumar 2019-6-18

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/467638-how-to-apply-the-boundary-conditons-to-the-mass-and-stiffness-matrices#answer_379706

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

Regards,

Ravi

3 个评论
显示 1更早的评论隐藏 1更早的评论

Ravi Kumar 2019-6-19

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

Sravan Kumar Putta 2021-2-25

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?

请先登录，再进行评论。