1. Resolução da Equação de Condução de Calor 3D por Elementos Finitos

Formulación Matemática:

\\(\\frac{\\partial u}{\\partial t} = \\alpha \\nabla^2 u + g(x,y,z,t)\\), onde \\(\\alpha\\) é o coeficiente de difusão térmica.

%% Programa FEM para Equação de Condução de Calor 3D
% Elementos tetraédricos lineares
clear; clc; close all;

%% Definição de Parâmetros
difusao_termica = 1.2;
tempo_final = 2.0;
passo_tempo = 0.005;
n_passos = ceil(tempo_final / passo_tempo);

%% Geração de Malha Tetraédrica (Domínio Cúbico)
xmin = -0.5; xmax = 0.5;
ymin = -0.5; ymax = 0.5;
zmin = -0.5; zmax = 0.5;
div_x = 6; div_y = 6; div_z = 6;

% Coordenadas dos nós
[Xg, Yg, Zg] = ndgrid(linspace(xmin, xmax, div_x+1), ...
                       linspace(ymin, ymax, div_y+1), ...
                       linspace(zmin, zmax, div_z+1));
coord_nos = [Xg(:), Yg(:), Zg(:)];
total_nos = size(coord_nos, 1);

% Construção dos elementos tetraédricos
lista_tetra = [];
for i_idx = 1:div_x
    for j_idx = 1:div_y
        for k_idx = 1:div_z
            % Índices dos vértices do hexaedro
            v1 = (i_idx-1)*(div_y+1)*(div_z+1) + (j_idx-1)*(div_z+1) + k_idx;
            v2 = v1 + 1;
            v3 = v1 + (div_z+1);
            v4 = v3 + 1;
            v5 = v1 + (div_y+1)*(div_z+1);
            v6 = v2 + (div_y+1)*(div_z+1);
            v7 = v3 + (div_y+1)*(div_z+1);
            v8 = v4 + (div_y+1)*(div_z+1);
            
            % Subdivisão em seis tetraedros
            lista_tetra = [lista_tetra; v1 v2 v3 v7];
            lista_tetra = [lista_tetra; v2 v3 v4 v7];
            lista_tetra = [lista_tetra; v1 v5 v6 v7];
            lista_tetra = [lista_tetra; v1 v3 v5 v7];
            lista_tetra = [lista_tetra; v3 v4 v7 v8];
            lista_tetra = [lista_tetra; v3 v6 v7 v8];
        end
    end
end
num_elementos = size(lista_tetra, 1);

%% Condição Inicial
solucao_inicial = zeros(total_nos, 1);
for n = 1:total_nos
    xn = coord_nos(n,1); yn = coord_nos(n,2); zn = coord_nos(n,3);
    solucao_inicial(n) = exp(-xn^2 - yn^2 - zn^2);
end

%% Identificação de Nós de Fronteira (Dirichlet)
tol_pos = 1e-8;
fronteira_idx = find(abs(coord_nos(:,1)-xmin) < tol_pos | abs(coord_nos(:,1)-xmax) < tol_pos | ...
                     abs(coord_nos(:,2)-ymin) < tol_pos | abs(coord_nos(:,2)-ymax) < tol_pos | ...
                     abs(coord_nos(:,3)-zmin) < tol_pos | abs(coord_nos(:,3)-zmax) < tol_pos);
valor_fronteira = zeros(size(fronteira_idx));

%% Montagem das Matrizes Globais
Matriz_rigidez = sparse(total_nos, total_nos);
Matriz_massa = sparse(total_nos, total_nos);
Vetor_forca = zeros(total_nos, 1);

for elem = 1:num_elementos
    nos_elem = lista_tetra(elem, :);
    coords_elem = coord_nos(nos_elem, :);
    
    [Ke_elem, Me_elem, Fe_elem] = calc_matrizes_tetraedro(coords_elem, difusao_termica);
    
    % Acoplamento na matriz global
    Matriz_rigidez(nos_elem, nos_elem) = Matriz_rigidez(nos_elem, nos_elem) + Ke_elem;
    Matriz_massa(nos_elem, nos_elem) = Matriz_massa(nos_elem, nos_elem) + Me_elem;
    Vetor_forca(nos_elem) = Vetor_forca(nos_elem) + Fe_elem;
end

%% Integração Temporal (Método θ, θ=0.5 para Crank-Nicolson)
param_theta = 0.6;
u_atual = solucao_inicial;

% Pré-cálculo das matrizes do sistema
Matriz_A = Matriz_massa + param_theta * passo_tempo * Matriz_rigidez;
Matriz_B = Matriz_massa - (1 - param_theta) * passo_tempo * Matriz_rigidez;

for passo = 1:n_passos
    instante = passo * passo_tempo;
    
    % Termo fonte dependente do tempo (exemplo: seno)
    F_temporal = Vetor_forca * sin(2*pi*instante);
    
    % Lado direito do sistema linear
    lado_dir = Matriz_B * u_atual + passo_tempo * ((1-param_theta)*F_temporal + param_theta*F_temporal);
    
    % Aplicar condições de contorno
    Matriz_corrigida = Matriz_A;
    rhs_corrigido = lado_dir;
    
    for idx = 1:length(fronteira_idx)
        no_fixo = fronteira_idx(idx);
        Matriz_corrigida(no_fixo, :) = 0;
        Matriz_corrigida(no_fixo, no_fixo) = 1;
        rhs_corrigido(no_fixo) = valor_fronteira(idx);
    end
    
    % Resolver sistema linear
    u_atual = Matriz_corrigida \ rhs_corrigido;
    
    if mod(passo, 20) == 0
        fprintf('Passo %d concluído (tempo %.4f)\n', passo, instante);
    end
end

%% Visualização
plano_z = 0.0;
indices_plano = find(abs(coord_nos(:,3) - plano_z) < 0.02);

figure;
scatter(coord_nos(indices_plano,1), coord_nos(indices_plano,2), 30, u_atual(indices_plano), 'filled');
colorbar;
title('Campo de temperatura no plano z=0');
xlabel('Coordenada X');
ylabel('Coordenada Y');
grid on;
axis equal;

%% Função de Cálculo das Matrizes Elementares
function [Ke, Me, Fe] = calc_matrizes_tetraedro(vertices, coef_dif)
    % vertices: matriz 4x3 com coordenadas dos nós do tetraedro
    
    % Volume do tetraedro
    V = abs(det([ones(4,1), vertices])) / 6;
    
    % Gradientes das funções de forma no espaço natural
    G_nat = [-1 1 0 0;
             -1 0 1 0;
             -1 0 0 1]';
    
    % Jacobiano e sua inversa
    J = G_nat' * vertices;
    inv_J = inv(J);
    
    % Gradientes no espaço físico
    G_fis = inv_J' * G_nat';
    
    % Matriz de rigidez elementar
    Ke = coef_dif * V * (G_fis * G_fis');
    
    % Matriz de massa elementar (fórmula fechada para linear)
    Me = (V/20) * (ones(4,4) + diag([1 1 1 1]));
    
    % Vetor de forças elementar (fonte unitária)
    Fe = (V/4) * ones(4,1);
end

2. Resolução da Equação de Poisson 3D por Elementos Finitos

Formulación Matemática:

\\(-\\nabla^2 u = f(x,y,z)\\), com condições de contorno adequadas.

%% Programa FEM para Equação de Poisson 3D
clear; clc; close all;

%% Definição do Problema
% Fonte f(x,y,z)
funcao_fonte = @(x,y,z) 6 * pi^2 * cos(pi*x) .* cos(pi*y) .* cos(pi*z);

%% Geração de Malha Regular
Lx = 1.0; Ly = 1.0; Lz = 1.0;
nx = 8; ny = 8; nz = 8;

% Nós da malha
[xn, yn, zn] = ndgrid(linspace(0, Lx, nx+1), linspace(0, Ly, ny+1), linspace(0, Lz, nz+1));
nos_poisson = [xn(:), yn(:), zn(:)];
total_nos_p = size(nos_poisson, 1);

% Elementos tetraédricos (subdivisão alternativa)
tetra_poisson = [];
for ii = 1:nx
    for jj = 1:ny
        for kk = 1:nz
            % Vértices do hexaedro local
            A = (ii-1)*(ny+1)*(nz+1) + (jj-1)*(nz+1) + kk;
            B = A + 1;
            C = A + (nz+1);
            D = C + 1;
            E = A + (ny+1)*(nz+1);
            F = B + (ny+1)*(nz+1);
            G = C + (ny+1)*(nz+1);
            H = D + (ny+1)*(nz+1);
            
            % Quatro tetraedros por hexaedro
            tetra_poisson = [tetra_poisson; A B D E];
            tetra_poisson = [tetra_poisson; B D E H];
            tetra_poisson = [tetra_poisson; A D E G];
            tetra_poisson = [tetra_poisson; D E G H];
        end
    end
end
num_elem_p = size(tetra_poisson, 1);

%% Condições de Contorno (Dirichlet Homogêneo)
tol_c = 1e-12;
contorno_idx = find(abs(nos_poisson(:,1)) < tol_c | abs(nos_poisson(:,1)-Lx) < tol_c | ...
                    abs(nos_poisson(:,2)) < tol_c | abs(nos_poisson(:,2)-Ly) < tol_c | ...
                    abs(nos_poisson(:,3)) < tol_c | abs(nos_poisson(:,3)-Lz) < tol_c);
valor_contorno = zeros(length(contorno_idx), 1);

%% Montagem do Sistema Global
K_global_p = sparse(total_nos_p, total_nos_p);
F_global_p = zeros(total_nos_p, 1);

for e = 1:num_elem_p
    nos_e = tetra_poisson(e, :);
    coords_e = nos_poisson(nos_e, :);
    
    [Ke_e, Fe_e] = elem_poisson_matrizes(coords_e, funcao_fonte);
    
    K_global_p(nos_e, nos_e) = K_global_p(nos_e, nos_e) + Ke_e;
    F_global_p(nos_e) = F_global_p(nos_e) + Fe_e;
end

%% Imposição das Condições de Contorno
K_fixa = K_global_p;
F_fixa = F_global_p;

for idx = 1:length(contorno_idx)
    no_cte = contorno_idx(idx);
    K_fixa(no_cte, :) = 0;
    K_fixa(no_cte, no_cte) = 1;
    F_fixa(no_cte) = valor_contorno(idx);
end

%% Resolução do Sistema Linear
solucao_p = K_fixa \ F_fixa;

%% Solução Analítica para Verificação
solucao_exata = cos(pi*nos_poisson(:,1)) .* cos(pi*nos_poisson(:,2)) .* cos(pi*nos_poisson(:,3));

%% Análise de Erro
erro_L2 = sqrt(mean((solucao_p - solucao_exata).^2));
erro_max = max(abs(solucao_p - solucao_exata));
fprintf('Erro L2: %e\n', erro_L2);
fprintf('Erro máximo: %e\n', erro_max);

%% Gráficos Comparativos
figure('Units','normalized','Position',[0.1 0.1 0.8 0.5]);

subplot(1,2,1);
scatter3(nos_poisson(:,1), nos_poisson(:,2), nos_poisson(:,3), 15, solucao_p, 'filled');
colorbar; title('Solução Numérica'); xlabel('X'); ylabel('Y'); zlabel('Z'); view(30,25);

subplot(1,2,2);
scatter3(nos_poisson(:,1), nos_poisson(:,2), nos_poisson(:,3), 15, abs(solucao_p-solucao_exata), 'filled');
colorbar; title('Distribuição do Erro'); xlabel('X'); ylabel('Y'); zlabel('Z'); view(30,25);

%% Função para Cálculo das Matrizes de Poisson
function [Ke, Fe] = elem_poisson_matrizes(vert, src_func)
    % Cálculo do volume
    V = abs(det([ones(4,1), vert])) / 6;
    
    % Gradientes naturais
    G_nat = [-1 1 0 0;
             -1 0 1 0;
             -1 0 0 1]';
    
    % Jacobiano
    J = G_nat' * vert;
    inv_J = inv(J);
    G_fis = inv_J' * G_nat';
    
    % Matriz de rigidez
    Ke = V * (G_fis * G_fis');
    
    % Vetor de forças (integração por ponto central)
    baricentro = mean(vert, 1);
    valor_f = src_func(baricentro(1), baricentro(2), baricentro(3));
    Fe = (V/4) * valor_f * ones(4,1);
end

3. Uso Avançado: Toolbox PDE do MATLAB

%% Exemplo com a Toolbox PDE do MATLAB
clear; clc; close all;

modelo_pde = createpde();

% Geometria de esfera sólida
g = multisphere(1);
modelo_pde.Geometry = g;

% Geração de malha adaptativa
gerarMalha(modelo_pde, 'Hmax', 0.15, 'Hmin', 0.05);

% Visualizar malha
pdeplot3D(modelo_pde);
title('Malha FEM 3D gerada');

% Especificação dos coeficientes
specifyCoefficients(modelo_pde, 'm', 0, 'd', 0, 'c', 1, 'a', 0, 'f', @(loc,state) loc.x.^2 + loc.y.^2);

% Condição de contorno (Neumann nula)
applyBoundaryCondition(modelo_pde, 'neumann', 'Face', 1:g.NumFaces, 'g', 0, 'q', 0);

% Resolver
res_pde = solvepde(modelo_pde);
u_pde = res_pde.NodalSolution;

figure;
pdeplot3D(modelo_pde, 'ColorMapData', u_pde);
title('Campo de Solução da EDP');
colorbar;

4. Refinamento Adaptativo de Malha

%% Algoritmo de Refinamento Adaptativo
function refinamento_adaptativo_3d()
    % Malha inicial grosseira
    [nos_ini, tetras_ini] = criar_malha_cubica(1);
    
    nivel_maximo = 5;
    tolerancia_ref = 0.01;
    
    for nivel = 1:nivel_maximo
        % Resolver no nível atual
        [u_nivel, estimativa_erro] = resolver_sistema(nos_ini, tetras_ini);
        
        % Marcar elementos para refinamento
        limiar = tolerancia_ref * max(estimativa_erro);
        marcar_refinamento = estimativa_erro > limiar;
        
        if ~any(marcar_refinamento)
            fprintf('Convergência atingida no nível %d\n', nivel);
            break;
        end
        
        % Refinar a malha
        [nos_ini, tetras_ini] = subdividir_elementos(nos_ini, tetras_ini, marcar_refinamento);
        fprintf('Nível %d: %d elementos refinados\n', nivel, sum(marcar_refinamento));
    end
    
    % Visualizar solução final
    figure;
    tetramesh(tetras_ini, nos_ini, u_nivel);
    title('Solução após Refinamento Adaptativo');
    colormap jet;
    colorbar;
end

Considerações Técnicas:

• A geração de malha 3D é crítica; para geometrias complexas, recomendam-se ferramentas como Gmsh ou distmesh.
• Elementos tetraédricos lineares oferecem bom compromisso entre precisão e custo computacional.
• O uso de matrizes esparsas é essencial para eficiência em problemas de grande porte.
• Para problemas transientes, métodos implícitos como Crenk-Nicolson proporcionam estabilidade.
• Em aplicações práticos, considere validação contra soluções analíticas ou dados experimentais.