fem_runtime.h

#ifndef MATRIX_FREE_MICROMOUSE_FEM_RUNTIME_H
#define MATRIX_FREE_MICROMOUSE_FEM_RUNTIME_H

typedef struct FemDeviceBuffers {
    CUdeviceptr rest;
    CUdeviceptr position;
    CUdeviceptr velocity;
    CUdeviceptr force;
    CUdeviceptr mass;
    CUdeviceptr fixed_mask;
    CUdeviceptr region;
    CUdeviceptr elements;
    CUdeviceptr element_material;
    CUdeviceptr materials;
    CUdeviceptr ties;
    CUdeviceptr tie_delta;
    CUdeviceptr gradients;
    CUdeviceptr dvol;
    CUdeviceptr element_metrics;
    CUdeviceptr history;
    CUdeviceptr point_dissipation;
} FemDeviceBuffers;

static void fem_device_free(CudaState *cuda, FemDeviceBuffers *device) {
#define FEM_FREE(field) do { if (device->field) cuda->cuMemFree(device->field); } while (0)
    FEM_FREE(point_dissipation);
    FEM_FREE(history);
    FEM_FREE(element_metrics);
    FEM_FREE(dvol);
    FEM_FREE(gradients);
    FEM_FREE(tie_delta);
    FEM_FREE(ties);
    FEM_FREE(materials);
    FEM_FREE(element_material);
    FEM_FREE(elements);
    FEM_FREE(region);
    FEM_FREE(fixed_mask);
    FEM_FREE(mass);
    FEM_FREE(force);
    FEM_FREE(velocity);
    FEM_FREE(position);
    FEM_FREE(rest);
#undef FEM_FREE
    memset(device, 0, sizeof(*device));
}

static CUresult fem_device_allocate(CudaState *cuda, FemDeviceBuffers *device,
                                    const FemModel *model) {
    CUresult result = 0;
    SIZE_T vector_bytes = (SIZE_T)model->node_count * 3u * sizeof(float);
#define FEM_ALLOC(field, bytes) \
    do { if ((result = cuda->cuMemAlloc(&device->field, (bytes)))) return result; } while (0)
    FEM_ALLOC(rest, vector_bytes);
    FEM_ALLOC(position, vector_bytes);
    FEM_ALLOC(velocity, vector_bytes);
    FEM_ALLOC(force, vector_bytes);
    FEM_ALLOC(mass, (SIZE_T)model->node_count * sizeof(float));
    FEM_ALLOC(fixed_mask, (SIZE_T)model->node_count * sizeof(u8));
    FEM_ALLOC(region, (SIZE_T)model->node_count * sizeof(u8));
    FEM_ALLOC(elements, (SIZE_T)model->element_count * 8u * sizeof(u32));
    FEM_ALLOC(element_material, (SIZE_T)model->element_count * sizeof(u32));
    FEM_ALLOC(materials, FEM_MAX_MATERIALS * FEM_MATERIAL_GPU_FLOATS * sizeof(float));
    FEM_ALLOC(ties, (SIZE_T)model->tie_count * 2u * sizeof(u32));
    FEM_ALLOC(tie_delta, (SIZE_T)model->tie_count * 3u * sizeof(float));
    FEM_ALLOC(gradients, (SIZE_T)model->element_count *
              FEM_GRADIENT_FLOATS_PER_ELEMENT * sizeof(float));
    FEM_ALLOC(dvol, (SIZE_T)model->element_count * 8u * sizeof(float));
    FEM_ALLOC(element_metrics, (SIZE_T)model->element_count *
              FEM_ELEMENT_METRIC_FLOATS * sizeof(float));
    FEM_ALLOC(history, (SIZE_T)model->element_count * FEM_GAUSS_POINTS *
              FEM_HISTORY_FLOATS_PER_POINT * sizeof(float));
    FEM_ALLOC(point_dissipation, (SIZE_T)model->element_count *
              FEM_GAUSS_POINTS * sizeof(float));
#undef FEM_ALLOC
    return 0;
}

static CUresult fem_device_upload(CudaState *cuda, FemDeviceBuffers *device,
                                  const FemModel *model) {
    CUresult result = 0;
    SIZE_T vector_bytes = (SIZE_T)model->node_count * 3u * sizeof(float);
    float materials[FEM_MAX_MATERIALS * FEM_MATERIAL_GPU_FLOATS];
    for (u32 material = 0u; material < FEM_MAX_MATERIALS; ++material) {
        u32 base = material * FEM_MATERIAL_GPU_FLOATS;
        materials[base] = model->materials[material].lambda;
        materials[base + 1u] = model->materials[material].mu_equilibrium;
        for (u32 branch = 0u; branch < FEM_MAXWELL_BRANCHES; ++branch) {
            materials[base + 2u + branch * 2u] =
                model->materials[material].branch_mu[branch];
            materials[base + 3u + branch * 2u] =
                model->materials[material].branch_tau[branch];
        }
    }
#define FEM_UPLOAD(field, source, bytes) \
    do { if ((result = cuda_copy_h2d(cuda, device->field, (source), (bytes)))) return result; } while (0)
    FEM_UPLOAD(rest, model->rest, vector_bytes);
    FEM_UPLOAD(position, model->position, vector_bytes);
    FEM_UPLOAD(velocity, model->velocity, vector_bytes);
    FEM_UPLOAD(force, model->force, vector_bytes);
    FEM_UPLOAD(mass, model->mass, (SIZE_T)model->node_count * sizeof(float));
    FEM_UPLOAD(fixed_mask, model->fixed_mask, (SIZE_T)model->node_count * sizeof(u8));
    FEM_UPLOAD(region, model->node_region, (SIZE_T)model->node_count * sizeof(u8));
    FEM_UPLOAD(elements, model->elements,
               (SIZE_T)model->element_count * 8u * sizeof(u32));
    FEM_UPLOAD(element_material, model->element_material,
               (SIZE_T)model->element_count * sizeof(u32));
    FEM_UPLOAD(materials, materials, sizeof(materials));
    FEM_UPLOAD(ties, model->ties, (SIZE_T)model->tie_count * 2u * sizeof(u32));
    FEM_UPLOAD(tie_delta, model->tie_delta,
               (SIZE_T)model->tie_count * 3u * sizeof(float));
    FEM_UPLOAD(gradients, model->gradients,
               (SIZE_T)model->element_count * FEM_GRADIENT_FLOATS_PER_ELEMENT *
               sizeof(float));
    FEM_UPLOAD(dvol, model->dvol,
               (SIZE_T)model->element_count * 8u * sizeof(float));
    FEM_UPLOAD(history, model->history, (SIZE_T)model->element_count *
               FEM_GAUSS_POINTS * FEM_HISTORY_FLOATS_PER_POINT * sizeof(float));
    FEM_UPLOAD(point_dissipation, model->point_dissipation,
               (SIZE_T)model->element_count * FEM_GAUSS_POINTS * sizeof(float));
#undef FEM_UPLOAD
    return 0;
}

static CUresult fem_launch_force(CudaState *cuda, const FemDeviceBuffers *device,
                                 u32 element_count, float dt,
                                 u32 update_history) {
    void *arguments[14] = {
        (void *)&device->rest, (void *)&device->position,
        (void *)&device->elements,
        (void *)&device->element_material, (void *)&device->materials,
        (void *)&device->gradients, (void *)&device->dvol,
        (void *)&device->force, (void *)&device->element_metrics,
        (void *)&device->history, (void *)&device->point_dissipation,
        &dt, &update_history, &element_count
    };
    u32 blocks = (element_count + 127u) / 128u;
    return cuda->cuLaunchKernel(cuda->fem_force, blocks, 1u, 1u,
                                128u, 1u, 1u, 0u, 0, arguments, 0);
}

static CUresult fem_launch_clear(CudaState *cuda, const FemDeviceBuffers *device,
                                 u32 node_count, const FemLoadCase *load,
                                 float load_scale) {
    float body_ax = load->body_ax, body_ay = load->body_ay;
    float body_az = load->body_az, floor_k = load->floor_stiffness;
    float floor_c = load->floor_damping, friction_mu = load->friction_mu;
    void *arguments[14] = {
        (void *)&device->rest, (void *)&device->position,
        (void *)&device->velocity, (void *)&device->mass,
        (void *)&device->region, (void *)&device->force, &node_count,
        &body_ax, &body_ay, &body_az, &floor_k, &floor_c, &friction_mu,
        &load_scale
    };
    u32 blocks = (node_count + 255u) / 256u;
    return cuda->cuLaunchKernel(cuda->fem_clear, blocks, 1u, 1u,
                                256u, 1u, 1u, 0u, 0, arguments, 0);
}

static CUresult fem_launch_ties(CudaState *cuda, const FemDeviceBuffers *device,
                                u32 tie_count, const FemLoadCase *load) {
    float tie_k = load->tie_stiffness, tie_c = load->tie_damping;
    void *arguments[8] = {
        (void *)&device->position, (void *)&device->velocity,
        (void *)&device->ties, (void *)&device->tie_delta,
        (void *)&device->force, &tie_count, &tie_k, &tie_c
    };
    u32 blocks = (tie_count + 127u) / 128u;
    return cuda->cuLaunchKernel(cuda->fem_tie, blocks, 1u, 1u,
                                128u, 1u, 1u, 0u, 0, arguments, 0);
}

static CUresult fem_launch_integrate(CudaState *cuda,
                                     const FemDeviceBuffers *device,
                                     u32 node_count, const FemLoadCase *load) {
    float dt = load->dt, damping = load->velocity_damping;
    void *arguments[9] = {
        (void *)&device->rest, (void *)&device->position,
        (void *)&device->velocity, (void *)&device->force,
        (void *)&device->mass, (void *)&device->fixed_mask,
        &node_count, &dt, &damping
    };
    u32 blocks = (node_count + 255u) / 256u;
    return cuda->cuLaunchKernel(cuda->fem_integrate, blocks, 1u, 1u,
                                256u, 1u, 1u, 0u, 0, arguments, 0);
}

static void fem_scale_geometry(FemModel *model, const SolverParams *params) {
    float track_scale = f_max((float)params->track_mm / 72.0f, 0.05f);
    float length_scale = f_max((float)params->chassis_length_mm / 90.0f, 0.05f);
    float wheel_scale = f_max((float)params->wheel_radius_mm / 12.0f, 0.05f);
    float half_track = 0.001f * 0.5f * (float)params->track_mm;
    float wheel_radius = 0.001f * (float)params->wheel_radius_mm;
    for (u32 node = 0u; node < model->node_count; ++node) {
        float x = model->rest[node * 3u + 0u];
        float y = model->rest[node * 3u + 1u];
        float z = model->rest[node * 3u + 2u];
        if (model->node_region[node] != FEM_REGION_LEFT_TIRE &&
            model->node_region[node] != FEM_REGION_RIGHT_TIRE) {
            x *= track_scale;
            y *= length_scale;
            z = wheel_radius * 1.5f + (z - 0.018f) * wheel_scale;
        } else {
            float center = model->node_region[node] == FEM_REGION_LEFT_TIRE ?
                           -half_track : half_track;
            float old_center = model->node_region[node] == FEM_REGION_LEFT_TIRE ?
                               -0.036f : 0.036f;
            x = center + (x - old_center) * wheel_scale;
            y *= wheel_scale;
            z = wheel_radius + (z - 0.012f) * wheel_scale;
        }
        model->rest[node * 3u + 0u] = x;
        model->rest[node * 3u + 1u] = y;
        model->rest[node * 3u + 2u] = z;
        model->position[node * 3u + 0u] = x;
        model->position[node * 3u + 1u] = y;
        model->position[node * 3u + 2u] = z;
    }
    for (u32 tie = 0u; tie < model->tie_count; ++tie) {
        u32 a = model->ties[tie * 2u + 0u];
        u32 b = model->ties[tie * 2u + 1u];
        for (u32 axis = 0u; axis < 3u; ++axis)
            model->tie_delta[tie * 3u + axis] =
                model->rest[a * 3u + axis] - model->rest[b * 3u + axis];
    }
}

static void fem_normalize_mass(FemModel *model, float target_mass) {
    float total = 0.0f;
    for (u32 node = 0u; node < model->node_count; ++node) total += model->mass[node];
    float scale = target_mass / f_max(total, 0.000001f);
    for (u32 node = 0u; node < model->node_count; ++node) model->mass[node] *= scale;
}

static void fem_route_load_envelope(const AppState *app, FemLoadCase *load,
                                    float *peak_longitudinal,
                                    float *peak_lateral) {
    float longitudinal = 0.0f;
    float lateral = 0.0f;
    float sample_dt = f_max((float)app->params.dt_us * 0.000004f, 0.000001f);
    for (u32 i = 0u; i < app->trajectory_count; ++i) {
        const TrajectorySample *sample = &app->trajectory[i];
        lateral = f_max(lateral,
                        f_abs(sample->speed_m_s * sample->yaw_rate_rad_s));
        if (i) {
            float acceleration = f_abs(sample->speed_m_s -
                                       app->trajectory[i - 1u].speed_m_s) / sample_dt;
            longitudinal = f_max(longitudinal, acceleration);
        }
    }
    longitudinal = f_clamp(longitudinal, 0.0f, 30.0f);
    lateral = f_clamp(lateral, 0.0f, 30.0f);
    load->body_ax = -lateral;
    load->body_ay = -longitudinal;
    load->body_az = -49.05f;
    *peak_longitudinal = longitudinal;
    *peak_lateral = lateral;
}

static float fem_effective_wheel_radius(const FemModel *model) {
    float total = 0.0f;
    u32 count = 0u;
    for (u32 tire = 0u; tire < 2u; ++tire) {
        u8 region = tire ? FEM_REGION_RIGHT_TIRE : FEM_REGION_LEFT_TIRE;
        u32 base = model->region_first[region];
        if (base == ~0u) continue;
        float center_y = 0.0f, center_z = 0.0f;
        u32 center_count = 0u;
        for (u32 ix = 0u; ix <= FEM_TIRE_NX; ++ix)
            for (u32 it = 0u; it < FEM_TIRE_NT; ++it) {
                u32 node = base + (ix * FEM_TIRE_NT + it) *
                                  (FEM_TIRE_NR + 1u);
                center_y += model->position[node * 3u + 1u];
                center_z += model->position[node * 3u + 2u];
                ++center_count;
            }
        center_y /= (float)center_count;
        center_z /= (float)center_count;
        for (u32 ix = 0u; ix <= FEM_TIRE_NX; ++ix)
            for (u32 it = 0u; it < FEM_TIRE_NT; ++it) {
                u32 node = base + (ix * FEM_TIRE_NT + it) *
                                  (FEM_TIRE_NR + 1u) + FEM_TIRE_NR;
                float dy = model->position[node * 3u + 1u] - center_y;
                float dz = model->position[node * 3u + 2u] - center_z;
                total += f_sqrt(dy * dy + dz * dz);
                ++count;
            }
    }
    return count ? total / (float)count : 0.0f;
}

static int run_fem_certification(AppState *app) {
    CudaState *cuda = &app->cuda;
    FemDeviceBuffers device;
    FemModel model;
    FemMetrics metrics;
    FemLoadCase load;
    CUresult result = 0;
    void *workspace = 0;
    memset(&device, 0, sizeof(device));
    memset(&g_fem_result, 0, sizeof(g_fem_result));
    g_fem_result.failure_reason = FEM_FAIL_WORKSPACE;
    workspace = VirtualAlloc(0, FEM_WORKSPACE_BYTES + 128u,
                             MEM_COMMIT | MEM_RESERVE, PAGE_READWRITE);
    if (!workspace || !fem_bind_workspace(&model, workspace,
                                           FEM_WORKSPACE_BYTES + 128u)) {
        wcopy(cuda->error, ARRAY_COUNT(cuda->error), L"全阶有限元主机工作区分配失败");
        goto fail_no_code;
    }
    if (!fem_build_mouse_model(&model)) {
        g_fem_result.failure_reason = FEM_FAIL_MESH_CAPACITY;
        wcopy(cuda->error, ARRAY_COUNT(cuda->error), L"全阶有限元网格容量不足");
        goto fail_no_code;
    }
    fem_scale_geometry(&model, &app->params);
    if (!fem_precompute_reference(&model)) {
        g_fem_result.failure_reason = FEM_FAIL_REFERENCE_JACOBIAN;
        wcopy(cuda->error, ARRAY_COUNT(cuda->error), L"全阶有限元参考构形 Jacobian 非正");
        goto fail_no_code;
    }
    if (!fem_validate_topology_and_materials(&model)) {
        g_fem_result.failure_reason = FEM_FAIL_TOPOLOGY;
        wcopy(cuda->error, ARRAY_COUNT(cuda->error),
              L"全车网格拓扑、质量或 Maxwell 分支审计失败");
        goto fail_no_code;
    }
    fem_normalize_mass(&model, 0.001f * (float)app->params.mass_g);
    fem_default_load_case(&load);
    fem_route_load_envelope(app, &load, &g_fem_result.peak_longitudinal_accel,
                            &g_fem_result.peak_lateral_accel);
    g_fem_result.node_count = model.node_count;
    g_fem_result.element_count = model.element_count;
    g_fem_result.dof_count = model.node_count * 3u;
    g_fem_result.active_dof_count = fem_active_dof_count(&model);
    g_fem_result.tie_count = model.tie_count;
    g_fem_result.gauss_point_count = model.element_count * 8u;
    g_fem_result.history_variable_count = model.element_count * 8u *
                                          FEM_HISTORY_FLOATS_PER_POINT;
    g_fem_result.material_count = FEM_MAX_MATERIALS;
    g_fem_result.steps = load.steps;

    {
        const float benchmark_deformation[9] = {
            1.08f, 0.04f, 0.00f,
            0.00f, 0.97f, 0.02f,
            0.00f, 0.00f, 0.95f
        };
        const float translation[3] = {0.001f, -0.002f, 0.003f};
        fem_apply_affine(&model, benchmark_deformation, translation);
    }
    memset(model.force, 0, (SIZE_T)model.node_count * 3u * sizeof(float));
    if (!fem_compute_internal_residual(&model, &metrics, 0, load.dt, 0)) {
        g_fem_result.failure_reason = FEM_FAIL_NONFINITE;
        wcopy(cuda->error, ARRAY_COUNT(cuda->error), L"CPU 有限应变参考残量出现非有限数");
        goto fail_no_code;
    }
    memcpy(model.reference_force, model.force,
           (SIZE_T)model.node_count * 3u * sizeof(float));

    if ((result = fem_device_allocate(cuda, &device, &model))) goto fail_cuda;
    if ((result = fem_device_upload(cuda, &device, &model))) goto fail_cuda;
    memset(model.force, 0, (SIZE_T)model.node_count * 3u * sizeof(float));
    if ((result = cuda_copy_h2d(cuda, device.force, model.force,
                                (SIZE_T)model.node_count * 3u * sizeof(float))))
        goto fail_cuda;
    if ((result = fem_launch_force(cuda, &device, model.element_count,
                                   load.dt, 0u))) goto fail_cuda;
    if ((result = cuda->cuCtxSynchronize())) goto fail_cuda;
    if ((result = cuda->cuMemcpyDtoH(model.force, device.force,
                (SIZE_T)model.node_count * 3u * sizeof(float)))) goto fail_cuda;
    if ((result = cuda->cuMemcpyDtoH(model.element_metrics, device.element_metrics,
                (SIZE_T)model.element_count * FEM_ELEMENT_METRIC_FLOATS *
                sizeof(float)))) goto fail_cuda;
    g_fem_result.force_relative_l2 = fem_force_relative_l2(
        model.reference_force, model.force, model.node_count * 3u);
    if (!(g_fem_result.force_relative_l2 == g_fem_result.force_relative_l2) ||
        g_fem_result.force_relative_l2 > 0.010f) {
        g_fem_result.failure_reason = FEM_FAIL_CPU_GPU_MISMATCH;
        wcopy(cuda->error, ARRAY_COUNT(cuda->error), L"CPU/GPU 有限元残量相对误差超过 1%");
        goto fail_no_code;
    }

    {
        u32 repetitions = 8u;
        u64 elapsed = 0u;
        do {
            u64 begin = GetTickCount64();
            for (u32 repeat = 0u; repeat < repetitions; ++repeat)
                if (!fem_compute_internal_residual(&model, &metrics, 0,
                                                    load.dt, 0)) {
                    g_fem_result.failure_reason = FEM_FAIL_NONFINITE;
                    wcopy(cuda->error, ARRAY_COUNT(cuda->error), L"CPU 有限元基准出现非有限数");
                    goto fail_no_code;
                }
            elapsed = GetTickCount64() - begin;
            if (elapsed < 160u && repetitions < 512u) repetitions *= 2u;
            else break;
        } while (1);
        g_fem_result.cpu_force_ms = (float)elapsed / (float)repetitions;
    }
    {
        u32 repetitions = 64u;
        u64 elapsed = 0u;
        do {
            u64 begin = GetTickCount64();
            for (u32 repeat = 0u; repeat < repetitions; ++repeat)
                if ((result = fem_launch_force(cuda, &device, model.element_count,
                                               load.dt, 0u)))
                    goto fail_cuda;
            if ((result = cuda->cuCtxSynchronize())) goto fail_cuda;
            elapsed = GetTickCount64() - begin;
            if (elapsed < 160u && repetitions < 4096u) repetitions *= 2u;
            else break;
        } while (1);
        g_fem_result.gpu_force_ms = (float)elapsed / (float)repetitions;
    }
    g_fem_result.speedup = g_fem_result.cpu_force_ms /
        f_max(g_fem_result.gpu_force_ms, 0.000001f);

    fem_reset_state(&model);
    if ((result = cuda_copy_h2d(cuda, device.position, model.position,
                (SIZE_T)model.node_count * 3u * sizeof(float)))) goto fail_cuda;
    if ((result = cuda_copy_h2d(cuda, device.velocity, model.velocity,
                (SIZE_T)model.node_count * 3u * sizeof(float)))) goto fail_cuda;
    if ((result = cuda_copy_h2d(cuda, device.history, model.history,
                (SIZE_T)model.element_count * FEM_GAUSS_POINTS *
                FEM_HISTORY_FLOATS_PER_POINT * sizeof(float)))) goto fail_cuda;
    if ((result = cuda_copy_h2d(cuda, device.point_dissipation,
                model.point_dissipation, (SIZE_T)model.element_count *
                FEM_GAUSS_POINTS * sizeof(float)))) goto fail_cuda;
    {
        u64 begin = GetTickCount64();
        for (u32 step = 0u; step < load.steps; ++step) {
            float scale = fem_load_scale(&load, step);
            if ((result = fem_launch_clear(cuda, &device, model.node_count,
                                           &load, scale))) goto fail_cuda;
            if ((result = fem_launch_ties(cuda, &device, model.tie_count,
                                          &load))) goto fail_cuda;
            if ((result = fem_launch_force(cuda, &device, model.element_count,
                                           load.dt, 1u)))
                goto fail_cuda;
            if ((result = fem_launch_integrate(cuda, &device, model.node_count,
                                               &load))) goto fail_cuda;
        }
        if ((result = cuda->cuCtxSynchronize())) goto fail_cuda;
        g_fem_result.gpu_transient_ms = (float)(GetTickCount64() - begin);
    }
    if ((result = cuda->cuMemcpyDtoH(model.position, device.position,
                (SIZE_T)model.node_count * 3u * sizeof(float)))) goto fail_cuda;
    if ((result = cuda->cuMemcpyDtoH(model.velocity, device.velocity,
                (SIZE_T)model.node_count * 3u * sizeof(float)))) goto fail_cuda;
    if ((result = cuda->cuMemcpyDtoH(model.force, device.force,
                (SIZE_T)model.node_count * 3u * sizeof(float)))) goto fail_cuda;
    if ((result = cuda->cuMemcpyDtoH(model.element_metrics, device.element_metrics,
                (SIZE_T)model.element_count * FEM_ELEMENT_METRIC_FLOATS *
                sizeof(float)))) goto fail_cuda;
    if ((result = cuda->cuMemcpyDtoH(model.point_dissipation,
                device.point_dissipation, (SIZE_T)model.element_count *
                FEM_GAUSS_POINTS * sizeof(float)))) goto fail_cuda;
    fem_reduce_element_metrics(&model, &metrics);
    g_fem_result.min_j = metrics.min_j;
    g_fem_result.max_j = metrics.max_j;
    g_fem_result.strain_energy = metrics.strain_energy;
    g_fem_result.max_von_mises_pa = metrics.max_von_mises;
    g_fem_result.dissipated_energy_j = metrics.dissipated_energy;
    g_fem_result.max_branch_stress_pa = metrics.max_branch_stress;
    g_fem_result.max_displacement_m = metrics.max_displacement;
    g_fem_result.max_lateral_displacement_m = 0.0f;
    for (u32 node = 0u; node < model.node_count; ++node)
        g_fem_result.max_lateral_displacement_m = f_max(
            g_fem_result.max_lateral_displacement_m,
            f_abs(model.position[node * 3u] - model.rest[node * 3u]));
    g_fem_result.effective_wheel_radius_m = fem_effective_wheel_radius(&model);
    if (!(metrics.min_j == metrics.min_j) || !(metrics.max_j == metrics.max_j) ||
        !(metrics.max_displacement == metrics.max_displacement)) {
        g_fem_result.failure_reason = FEM_FAIL_NONFINITE;
        wcopy(cuda->error, ARRAY_COUNT(cuda->error), L"GPU 全阶瞬态出现非有限数");
        goto fail_no_code;
    }
    if (metrics.min_j < 0.20f) {
        g_fem_result.failure_reason = FEM_FAIL_ELEMENT_INVERSION;
        wcopy(cuda->error, ARRAY_COUNT(cuda->error), L"GPU 全阶瞬态检测到单元翻转（J < 0.2）");
        goto fail_no_code;
    }
    if (metrics.max_displacement > 0.050f) {
        g_fem_result.failure_reason = FEM_FAIL_DISPLACEMENT_LIMIT;
        wcopy(cuda->error, ARRAY_COUNT(cuda->error), L"GPU 全阶瞬态位移超过 50 mm 安全门");
        goto fail_no_code;
    }
    g_fem_result.valid = 1u;
    g_fem_result.failure_reason = FEM_FAIL_NONE;
    fem_device_free(cuda, &device);
    VirtualFree(workspace, 0, MEM_RELEASE);
    return 1;

fail_cuda:
    g_fem_result.failure_reason = FEM_FAIL_CUDA;
    cuda_set_error(cuda, L"CUDA 全阶有限元阶段失败", result);
fail_no_code:
    fem_device_free(cuda, &device);
    if (workspace) VirtualFree(workspace, 0, MEM_RELEASE);
    return 0;
}

#endif