地盤技術者が自作FEMプログラム用にGUIを開発する話（第3回）

はじめに

第3回ではGUIに変形図、コンター図を表示させます。
対象の問題としては以下の図のように軟弱な土の上に構造物を載せる問題を考えます。
構造物の重さで地盤沈下が起きます。

メッシュ図表示

メッシュ図は以下のように表示させます。

ui.label(format!("nodes={} elements={}", nodes.len(), elements.len()));
    let available_size = ui.available_size();
    let (response, painter) = ui.allocate_painter(available_size, egui::Sense::hover());
    painter.rect_filled(response.rect, 0.0, egui::Color32::WHITE);

    if nodes.is_empty() || elements.is_empty() {
        return;
    }

    let rect = response.rect.shrink(40.0);
    let Some(transform) = MeshTransform::new(nodes, rect) else {
        return;
    };
    let stroke = egui::Stroke::new(1.5, egui::Color32::BLACK);

    for element in elements {
        let Some(p1) = node_position(nodes, element.n1, &transform) else {
            continue;
        };
        let Some(p2) = node_position(nodes, element.n2, &transform) else {
            continue;
        };
        let Some(p3) = node_position(nodes, element.n3, &transform) else {
            continue;
        };
        let Some(p4) = node_position(nodes, element.n4, &transform) else {
            continue;
        };
        painter.line_segment([p1, p2], stroke);
        painter.line_segment([p2, p3], stroke);
        painter.line_segment([p3, p4], stroke);
        painter.line_segment([p4, p1], stroke);

        if show_element_ids {
            let center = egui::pos2(
                (p1.x + p2.x + p3.x + p4.x) / 4.0,
                (p1.y + p2.y + p3.y + p4.y) / 4.0,
            );

            painter.text(
                center,
                egui::Align2::CENTER_CENTER,
                element.element_id.to_string(),
                egui::FontId::proportional(12.0),
                egui::Color32::DARK_GRAY,
            );
        }

        if show_node_ids {
            draw_node_id(&painter, p1, element.n1);
            draw_node_id(&painter, p2, element.n2);
            draw_node_id(&painter, p3, element.n3);
            draw_node_id(&painter, p4, element.n4);
        }
    }

node_positionはこんな感じです。

nodes
        .iter()
        .find(|node| node.node_id == node_id)
        .map(|node| transform.world_to_screen(node.x, node.y))

draw_node_idはこのように書きます。

let offset = egui::vec2(6.0, -6.0);
    painter.text(
        position + offset,
        egui::Align2::LEFT_BOTTOM,
        node_id.to_string(),
        egui::FontId::proportional(11.0),
        egui::Color32::from_rgb(100, 149, 237),
    );

FEMでは節点を作るので、その節点どうしを線で結ぶイメージです。
結果はこのようになります。

変形図作成

変形図は同じ用にメッシュ作成をしますが、初期節点ではなく変形後の節点をつなぐ必要があります。
今、解析結果（節点変位）を以下のような構造体に入れます。

pub struct NodeResult {
    pub step: usize,
    pub time: f64,
    pub node_id: usize,
    pub ux: f64,
    pub uy: f64,
}

後はこれを読み込んで先ほどと同じ用に出力させます。

ui.label(format!("nodes={} elements={}", nodes.len(), elements.len()));
    let available_size = ui.available_size();
    let (response, painter) = ui.allocate_painter(available_size, egui::Sense::hover());
    painter.rect_filled(response.rect, 0.0, egui::Color32::WHITE);

    if nodes.is_empty() || elements.is_empty() {
        return;
    }

    let rect = response.rect.shrink(40.0);
    let Some(transform) = MeshTransform::new(nodes, rect) else {
        return;
    };

    let initial_stroke = egui::Stroke::new(1.0, egui::Color32::LIGHT_GRAY);
    let deformed_stroke = egui::Stroke::new(1.5, egui::Color32::RED);

    for element in elements {
        let Some(p1) = node_position(nodes, element.n1, &transform) else {
            continue;
        };
        let Some(p2) = node_position(nodes, element.n2, &transform) else {
            continue;
        };
        let Some(p3) = node_position(nodes, element.n3, &transform) else {
            continue;
        };
        let Some(p4) = node_position(nodes, element.n4, &transform) else {
            continue;
        };
        painter.line_segment([p1, p2], initial_stroke);
        painter.line_segment([p2, p3], initial_stroke);
        painter.line_segment([p3, p4], initial_stroke);
        painter.line_segment([p4, p1], initial_stroke);
        let Some(dp1) = deformed_node_position_by_id(
            nodes,
            node_results,
            element.n1,
            deformation_scale,
            &transform,
            selected_step,
        ) else {
            continue;
        };
        let Some(dp2) = deformed_node_position_by_id(
            nodes,
            node_results,
            element.n2,
            deformation_scale,
            &transform,
            selected_step,
        ) else {
            continue;
        };
        let Some(dp3) = deformed_node_position_by_id(
            nodes,
            node_results,
            element.n3,
            deformation_scale,
            &transform,
            selected_step,
        ) else {
            continue;
        };
        let Some(dp4) = deformed_node_position_by_id(
            nodes,
            node_results,
            element.n4,
            deformation_scale,
            &transform,
            selected_step,
        ) else {
            continue;
        };
        painter.line_segment([dp1, dp2], deformed_stroke);
        painter.line_segment([dp2, dp3], deformed_stroke);
        painter.line_segment([dp3, dp4], deformed_stroke);
        painter.line_segment([dp4, dp1], deformed_stroke);

        if show_element_ids {
            let center = egui::pos2(
                (dp1.x + dp2.x + dp3.x + dp4.x) / 4.0,
                (dp1.y + dp2.y + dp3.y + dp4.y) / 4.0,
            );

            painter.text(
                center,
                egui::Align2::CENTER_CENTER,
                element.element_id.to_string(),
                egui::FontId::proportional(12.0),
                egui::Color32::DARK_GRAY,
            );
        }

        if show_node_ids {
            draw_node_id(&painter, dp1, element.n1);
            draw_node_id(&painter, dp2, element.n2);
            draw_node_id(&painter, dp3, element.n3);
            draw_node_id(&painter, dp4, element.n4);
        }
    }

deformed_node_position_by_idはnode_positionとほぼ同じです。
これで以下のような変形図が書けます。

コンター図作成

コンター図は最初に作ったメッシュの図形の中を色塗りするイメージです。
以下のように実装します。

ui.label(format!("nodes={} elements={}", nodes.len(), elements.len()));
    let available_size = ui.available_size();
    let (response, painter) = ui.allocate_painter(available_size, egui::Sense::hover());
    painter.rect_filled(response.rect, 0.0, egui::Color32::WHITE);

    if nodes.is_empty() || elements.is_empty() {
        return;
    }

    let rect = response.rect.shrink(40.0);
    let Some(transform) = MeshTransform::new(nodes, rect) else {
        return;
    };
    let values = collect_contour_values(
        contour_type,
        nodes,
        elements,
        node_results,
        element_results,
        internal_variable_results,
        None,
    );

    if values.is_empty() {
        return;
    }

    let vmin = values.iter().copied().fold(f64::INFINITY, f64::min);
    let vmax = values.iter().copied().fold(f64::NEG_INFINITY, f64::max);
    ui.label(format!(
        "{} min={}, max={}",
        contour_type.label(),
        contour::format_value(vmin),
        contour::format_value(vmax),
    ));
    let mesh_stroke = egui::Stroke::new(1.0, egui::Color32::BLACK);

    for element in elements {
        let Some(p1) = contour_node_position(
            nodes,
            node_results,
            element.n1,
            &transform,
            selected_step,
            show_deformed_contour,
            deformation_scale,
        ) else {
            continue;
        };
        let Some(p2) = contour_node_position(
            nodes,
            node_results,
            element.n2,
            &transform,
            selected_step,
            show_deformed_contour,
            deformation_scale,
        ) else {
            continue;
        };
        let Some(p3) = contour_node_position(
            nodes,
            node_results,
            element.n3,
            &transform,
            selected_step,
            show_deformed_contour,
            deformation_scale,
        ) else {
            continue;
        };
        let Some(p4) = contour_node_position(
            nodes,
            node_results,
            element.n4,
            &transform,
            selected_step,
            show_deformed_contour,
            deformation_scale,
        ) else {
            continue;
        };
        let Some(value) = contour::contour_element_value(
            contour_type,
            node_results,
            element_results,
            internal_variable_results,
            element,
            selected_step,
        ) else {
            continue;
        };
        let fill_color = contour::contour_color(value, vmin, vmax);
        painter.add(egui::Shape::convex_polygon(
            vec![p1, p2, p3, p4],
            fill_color,
            mesh_stroke,
        ));
    }

contour::contour_element_valueはこんな感じです。

match contour_type.source() {
        ContourSource::Node => {
            let v1 = node_contour_value(node_results, element.n1, contour_type, selected_step)?;
            let v2 = node_contour_value(node_results, element.n2, contour_type, selected_step)?;
            let v3 = node_contour_value(node_results, element.n3, contour_type, selected_step)?;
            let v4 = node_contour_value(node_results, element.n4, contour_type, selected_step)?;
            Some((v1 + v2 + v3 + v4) / 4.0)
        }
        ContourSource::Element => {
            selected_element_result(element_results, element.element_id, selected_step)
                .map(|result| element_contour_value(contour_type, result))
        }
        ContourSource::InternalVariable => selected_internal_variable_result(
            internal_variable_results,
            element.element_id,
            selected_step,
        )
        .map(|result| internal_variable_contour_value(contour_type, result)),
    }

結果はこんな感じです。

おわりに

コンター図まで作成できてだいぶポストプロセッサーらしくなりました。