Avaliação vetorizada das formas modais, zone frequências de excitação

attenuamas.py

# attenuamas.py
import numpy as np

class FoldingModel:
    # Características:
    # (1) Integração semi-implícita estável para as amplitudes modais
    # (2) Avaliação vetorizada das formas modais
    # (3) Cálculo de deslocamento e atenuação por elemento
    # (4) Métodos de teste unitário integrados para verificar consistência e estabilidade
    def __init__(self,
                 zone_center_z=0.0,
                 zone_radius=0.35,
                 base_fold_amp=0.02,
                 modes=None,
                 modal_stiffness=50.0,
                 modal_damping=8.0,
                 modal_mass=1.0,
                 max_local_strain=0.25,
                 strain_atten_sensitivity=2.0,
                 radial_scale=0.12,
                 max_amp_clip=0.2):
        # Parâmetros geométricos e físicos:
        self.zone_center_z = float(zone_center_z)   # centro em z da zona de influência
        self.zone_radius = float(zone_radius)       # raio característico da zona
        self.base_fold_amp = float(base_fold_amp)   # amplitude base de forçamento (m)
        self.modes = list(modes) if modes is not None else [(1, 0), (2, 0)]
        # Parâmetros do modelo modal (massa-mola-amortecimento):
        self.K = float(modal_stiffness)
        self.C = float(modal_damping)
        self.M = float(modal_mass)
        # Parâmetros de atenuação/strain:
        self.max_local_strain = float(max_local_strain)
        self.strain_atten_sensitivity = float(strain_atten_sensitivity)
        # Parâmetros numéricos:
        self.radial_scale = float(radial_scale)
        self.max_amp_clip = float(max_amp_clip)
        # Estado modal; cada modo tem [amplitude, velocidade]:
        self.n_modes = len(self.modes)
        self.modal_state = np.zeros((self.n_modes, 2), dtype=np.float64)

    def _modal_shape(self, elem_pos):
        # Calcula as formas modais para todos os modos e elementos
        # Retorna um array (n_modes, N) contendo valores adimensionais da forma
        
        # Estratégia:
        # (1) Componente azimutal: cos(m * phi)
        # (2) Envelope radial: gaussiana centrada em r0 multiplicada por potência r^n_rad
        # (3) Localização em z: gaussiana centrada em zone_center_z
        x = elem_pos[:, 0].astype(np.float64)
        y = elem_pos[:, 1].astype(np.float64)
        z = elem_pos[:, 2].astype(np.float64)
        r = np.sqrt(x * x + y * y)
        phi = np.arctan2(y, x)
        r0 = self.zone_radius
        radial_base = np.exp(-((r - r0) ** 2) / (2 * (self.radial_scale ** 2)))
        shapes = np.empty((self.n_modes, elem_pos.shape[0]), dtype=np.float64)
        for i, (m_az, n_rad) in enumerate(self.modes):
            # Componente azimutal: modo real com cosseno
            az = np.cos(m_az * phi)
            # Envelope radial, cee potência segura seguida de gaussiana
            if n_rad <= 0:
                radial = radial_base
            else:
                radial = (r / (r0 + 1e-12)) ** n_rad
                radial *= radial_base
            # Localização em z (gaussiana para limitar influência)
            zenv = np.exp(-((z - self.zone_center_z) ** 2) / (2 * (0.15 ** 2)))
            shapes[i, :] = az * radial * zenv
        return shapes

    def _compute_local_strain(self, elem_pos, displacement):
        # Estima uma proxy de strain local com base no deslocamento radial
        # Fórmula usada: strain ≈ |Δr| / r_local
        # Para cada elemento; projeta deslocamento na direção radial e normaliza por r
        x = elem_pos[:, 0]; y = elem_pos[:, 1]
        r = np.sqrt(x * x + y * y) + 1e-12
        rx = x / r; ry = y / r
        dr = displacement[:, 0] * rx + displacement[:, 1] * ry
        strain = np.abs(dr) / r
        return strain

    def evaluate(self, elem_pos, mach=8.0, t=0.0, dt_physical=None):
        # Calcula deslocamentos (N,3) e atenuações (N,) para as posições elem_pos

        # Parâmetros:
        # elem_pos: array (N,3) com coordenadas em metros
        # mach: fator cinemático controlador (escalado internamente)
        # t: tempo atual (s), usado para fase de forçamento
        # dt_physical: passo de tempo para integração modal (s); adrão: 1e-3 s

        # Retorna:
        # ~displacement: array (N,3) em metros
        # ~attenuation: array (N,) valores em [0.01, 1.0]
        N = elem_pos.shape[0]
        dt = 1e-3 if dt_physical is None else float(dt_physical)
        # Mapeamento estável do parâmetro mach para fator de forçamento
        forcing_factor = np.clip(mach / 10.0, 0.0, 2.0)
        # Cálculo vetorizado das formas modais
        shapes = self._modal_shape(elem_pos)  # (n_modes, N)
        # Frequência natural modal aproximada
        wn = np.sqrt(self.K / (self.M + 1e-18))
        # Passo modal semi-implícito (vetorizado)
        # ~proj_strength = norma L2 das formas por modo
        proj_strength = np.sqrt(np.mean(shapes * shapes, axis=1) + 1e-20)  # (n_modes,)
        # Frequências de excitação por modo (fracionárias de wn)
        if self.n_modes > 1:
            mode_idx = np.arange(self.n_modes)
            drive_freq = wn * (0.5 + 0.5 * (mode_idx / (self.n_modes - 1)))
        else:
            drive_freq = np.array([wn * 0.5])
        F_t = self.base_fold_amp * forcing_factor * proj_strength * np.sin(2.0 * np.pi * drive_freq * t)
        vel = self.modal_state[:, 1]
        amp = self.modal_state[:, 0]
        denom = self.M / dt + self.C
        rhs = self.M / dt * vel + F_t - self.K * amp
        vel_new = rhs / denom
        amp_new = amp + vel_new * dt
        amp_new = np.clip(amp_new, -self.max_amp_clip, self.max_amp_clip)
        # Salva estado modal atualizado
        self.modal_state[:, 0] = amp_new
        self.modal_state[:, 1] = vel_new
        # Monta deslocamento total por elemento
        x = elem_pos[:, 0].astype(np.float64)
        y = elem_pos[:, 1].astype(np.float64)
        z = elem_pos[:, 2].astype(np.float64)
        r = np.sqrt(x * x + y * y) + 1e-12
        rx = x / r; ry = y / r
        radial_contrib = (amp_new[:, None] * shapes)  # (n_modes, N)
        radial_total = np.sum(radial_contrib, axis=0)  # (N,)
        z_total = 0.02 * radial_total
        displ_total = np.empty((N, 3), dtype=np.float64)
        displ_total[:, 0] = radial_total * rx
        displ_total[:, 1] = radial_total * ry
        displ_total[:, 2] = z_total
        # Controle de magnitude global para evitar deformações excessivas
        norms = np.linalg.norm(displ_total, axis=1)
        max_disp = 0.25
        too_big = norms > max_disp
        if np.any(too_big):
            displ_total[too_big] *= (max_disp / norms[too_big])[:, None]
        # Suavização espacial via pesos gaussianos
        radial_weight = np.exp(-((r - self.zone_radius) ** 2) / (2 * (self.radial_scale ** 2)))
        zweight = np.exp(-((z - self.zone_center_z) ** 2) / (2 * (0.15 ** 2)))
        combined_weight = radial_weight * zweight
        displ_total *= combined_weight[:, None]
        # Restrição radial: elementos não podem penetrar abaixo de min_r
        r0 = self.zone_radius
        min_r = 0.6 * r0
        radial_disp_proj = displ_total[:, 0] * rx + displ_total[:, 1] * ry
        new_r = r + radial_disp_proj
        collide = new_r < min_r
        if np.any(collide):
            rx_c = rx[collide]; ry_c = ry[collide]
            excess = (min_r - new_r[collide])
            displ_total[collide, 0] += excess * rx_c
            displ_total[collide, 1] += excess * ry_c
        # Cálculo de atenuação cee combinando strain e efeito z-plane:
        strain = self._compute_local_strain(elem_pos, displ_total)
        s_norm = strain / (self.max_local_strain + 1e-12)
        atten = 1.0 / (1.0 + self.strain_atten_sensitivity * (s_norm ** 2))
        zdist = np.abs(elem_pos[:, 2] - self.zone_center_z)
        z_atten = 1.0 - 0.5 * np.exp(- (zdist ** 2) / (2 * (0.12 ** 2))) * combined_weight
        atten *= z_atten
        atten = np.clip(atten, 0.01, 1.0)
        return displ_total, atten

    # Métodos de teste embutidos:
    # ____________________________
    def run_basic_consistency_test(self):
        # Teste básico de consistência:
        
        # (1) Verifica formas de saída (N,3) e (N,)
        # (2) Verifica limites de deslocamento e atenuação
        # (3) Verifica reprodutibilidade: duas chamadas consecutivas sem avanço de tempo
        # (mesmo t, mesmo estado inicial) não devem produzir divergência abrupta
        
        # **Retorna True se passar, lança AssertionError caso contrário
        #   cria um conjunto simples de elementos: arranjo circular
        theta = np.linspace(0, 2*np.pi, 64, endpoint=False)
        r = np.full_like(theta, self.zone_radius)
        x = r * np.cos(theta); y = r * np.sin(theta); z = np.zeros_like(theta)
        elem_pos = np.stack([x, y, z], axis=1)
        # Salva estado modal para comparar reprodutibilidade
        saved_state = self.modal_state.copy()
        displ1, atten1 = self.evaluate(elem_pos, mach=8.0, t=0.123, dt_physical=1e-3)
        # Chamar novamente com mesmo t e verificar comportamento determinístico
        displ2, atten2 = self.evaluate(elem_pos, mach=8.0, t=0.123, dt_physical=1e-3)
        assert displ1.shape == (elem_pos.shape[0], 3), "Displacement tem shape inválido"
        assert atten1.shape == (elem_pos.shape[0],), "Attenuation tem shape inválido"
        # Limites de deslocamento e atenuação
        maxnorm = np.max(np.linalg.norm(displ1, axis=1))
        assert maxnorm <= 0.26, f"Deslocamento muito grande: {maxnorm}"
        assert np.all((atten1 >= 0.01) & (atten1 <= 1.0)), "Atenuação fora de [0.01,1]"
        # Reprodutibilidade: estados modais devem ter sido atualizados, mas sem NaNs
        assert np.all(np.isfinite(self.modal_state)), "Estado modal contém NaN/Inf"
        # Restaurar estado e verificar que chamada inicial era determinística
        self.modal_state = saved_state.copy()
        displ_restore, atten_restore = self.evaluate(elem_pos, mach=8.0, t=0.123, dt_physical=1e-3)
        # Comparar duas primeiras avaliações com mesmo estado inicial
        assert np.allclose(displ1, displ_restore, atol=1e-12) or np.allclose(atten1, atten_restore, atol=1e-12) \
            , "Comportamento não determinístico entre chamadas idênticas (verifique estado interno)"
        return True

    def run_stability_test(self, steps=200, dt=1e-3):
        # Teste de estabilidade temporal:
        # (1) Integra o sistema por 'steps' passos com forçamento constante
        # (2) Verifica que as amplitudes modais permanecem dentro do limite max_amp_clip
        #     e não explodem numericamente
        Retorna True se passar, AssertionError caso contrário.
        # ⌚Gera elementos aleatórios numa coroa para excitar modos:
        rng = np.random.default_rng(42)
        N = 512
        angles = rng.random(N) * 2*np.pi
        radii = self.zone_radius * (0.8 + 0.4 * rng.random(N))
        x = radii * np.cos(angles)
        y = radii * np.sin(angles)
        z = 0.02 * (rng.random(N) - 0.5)
        elem_pos = np.stack([x, y, z], axis=1)
        # 🗿Reset do estado modal para condições iniciais estáveis
        self.modal_state[:] = 0.0
        for k in range(steps):
            t = k * dt
            displ, atten = self.evaluate(elem_pos, mach=10.0, t=t, dt_physical=dt)
            amps = self.modal_state[:, 0]
            # Verificar limites numéricos e boundedness
            assert np.all(np.isfinite(amps)), "Estado modal contém NaN/Inf durante integração"
            assert np.max(np.abs(amps)) <= (self.max_amp_clip + 1e-6), \
                f"Amplitudes excederam max_amp_clip: {np.max(np.abs(amps))}"
            # Checar que atenuação permanece no intervalo
            assert np.all((atten >= 0.01) & (atten <= 1.0)), "Atenuação saiu do intervalo durante integração"
        return True

🗣🗿