'use client'
import React, { useRef, useMemo } from 'react';
import { Canvas, useFrame } from '@react-three/fiber';
import { OrbitControls, Sky, Line, Text } from '@react-three/drei';
import { useControls } from 'leva';
import * as THREE from 'three';

function NeutralAxisVisualization() {
  const meshRef = useRef();

  const { naRotation, naPosition, maxStrain, wireframe, showNeutralAxis, showStrainDiagram } = useControls({
    naRotation: { value: 0, min: -Math.PI, max: Math.PI, step: 0.01 },
    naPosition: { value: 0, min: -4, max: 4, step: 0.1 },
    maxStrain: { value: -0.0035 },
    wireframe: false,
    showNeutralAxis: true,
    showStrainDiagram: true
  });
  const { compression, tension } = getExtremeFibers(naRotation, naPosition, 2);
  const { distance: c, intersection } = getDistanceAndClosestPoint(compression.x, compression.y, naRotation, naPosition);
  console.log({ compression, intersection, c });

  const material = useMemo(() => {
    return new THREE.ShaderMaterial({
      vertexShader: `
        uniform float uNeutralAxisRotation;
        uniform float uNeutralAxisPosition;
        uniform float uC;
        uniform float uMaxStrain;
        varying vec3 vPosition;
        varying vec2 vUv;
        varying float vDistanceFromNA;
        varying float vStrain;
        
        void main() {
          vec3 newPosition = position;
          
          // Calculate distance from neutral axis
          // Rotate coordinate system by neutral axis rotation
          float cosAngle = cos(-1.0*uNeutralAxisRotation);
          float sinAngle = sin(-1.0*uNeutralAxisRotation);
          
          // Transform position to neutral axis coordinate system
          vec2 rotatedPos = vec2(
            position.x * cosAngle - position.y * sinAngle,
            position.x * sinAngle + position.y * cosAngle
          );
          
          // Distance from neutral axis (perpendicular distance)
          float distanceFromNA = rotatedPos.y - uNeutralAxisPosition;
          
          // Calculate strain using linear strain distribution: ε = (distance/c) * εmax
          float strain = (distanceFromNA / uC) * uMaxStrain;
        
          gl_Position = projectionMatrix * modelViewMatrix * vec4(newPosition, 1.0);
          vPosition = position;
          vUv = uv;
          vDistanceFromNA = distanceFromNA;
          vStrain = strain;
        }
      `,
      fragmentShader: `
        varying float vStrain;
      
        // Vibrant color spectrum function (blue -> cyan -> green -> yellow -> red)
        vec3 getStrainColor(float strain) {
          // Normalize strain to 0-1 range for color mapping
          float normalizedStrain = (strain + 1.0) * 0.5; // Map from [-1,1] to [0,1]
          
          vec3 color;
          
          if (normalizedStrain < 0.25) {
            // Deep blue to cyan
            float t = normalizedStrain * 4.0;
            color = mix(vec3(0.0, 0.0, 1.0), vec3(0.0, 1.0, 1.0), t);
          } else if (normalizedStrain < 0.5) {
            // Cyan to green
            float t = (normalizedStrain - 0.25) * 4.0;
            color = mix(vec3(0.0, 1.0, 1.0), vec3(0.0, 1.0, 0.0), t);
          } else if (normalizedStrain < 0.75) {
            // Green to yellow
            float t = (normalizedStrain - 0.5) * 4.0;
            color = mix(vec3(0.0, 1.0, 0.0), vec3(1.0, 1.0, 0.0), t);
          } else {
            // Yellow to red
            float t = (normalizedStrain - 0.75) * 4.0;
            color = mix(vec3(1.0, 1.0, 0.0), vec3(1.0, 0.0, 0.0), t);
          }
          
          return color;
        }
      
        void main() {
          // Normalize strain for color mapping (-maxStrain to +maxStrain -> -1 to +1)
          float normalizedStrain = clamp(vStrain / 0.003, -1.0, 1.0);
          vec3 color = getStrainColor(normalizedStrain);
          
          gl_FragColor = vec4(color, 1.0);
        }
      `,
      uniforms: {
        uNeutralAxisRotation: { value: naRotation },
        uNeutralAxisPosition: { value: naPosition },
        uC: { value: c },
        uMaxStrain: { value: maxStrain }
      },
      side: THREE.DoubleSide,
      wireframe,
    });
  }, [naRotation, naPosition, c, maxStrain, wireframe, showNeutralAxis]);

  useFrame((state) => {
    // Update uniforms only - no rotation
    if (material.uniforms) {
      material.uniforms.uNeutralAxisRotation.value = naRotation;
      material.uniforms.uNeutralAxisPosition.value = naPosition;
      material.uniforms.uC.value = c;
      material.uniforms.uMaxStrain.value = maxStrain;
    }
  });

  // Calculate neutral axis line endpoints
  const neutralAxisLine = useMemo(() => {
    const cosAngle = Math.cos(naRotation);
    const sinAngle = Math.sin(naRotation);

    // Line length extending across the geometry
    const lineLength = 4;

    // Calculate perpendicular offset for the neutral axis position
    const offsetX = -sinAngle * naPosition;
    const offsetY = cosAngle * naPosition;

    // Line endpoints
    const start = [
      offsetX - cosAngle * lineLength,
      offsetY - sinAngle * lineLength,
      0.01
    ];
    const end = [
      offsetX + cosAngle * lineLength,
      offsetY + sinAngle * lineLength,
      0.01
    ];

    return [start, end];
  }, [naRotation, naPosition]);

  // Calculate strain diagram line (perpendicular to NA, bounded by section geometry)
  const strainDiagramLine = useMemo(() => {
    const cosAngle = Math.cos(naRotation);
    const sinAngle = Math.sin(naRotation);

    // Perpendicular direction (rotated 90 degrees)
    const perpCos = -sinAngle;
    const perpSin = cosAngle;

    // Offset for neutral axis position
    const offsetX = -sinAngle * naPosition;
    const offsetY = cosAngle * naPosition;

    // Section bounds (assuming 4x4 geometry from -2 to +2)
    const sectionHalfWidth = 2.0;
    const sectionHalfHeight = 2.0;

    // Find intersection points with section boundaries
    const intersections = [];

    // Check intersection with each boundary
    const boundaries = [
      { normal: [1, 0], point: [sectionHalfWidth, 0] },    // right
      { normal: [-1, 0], point: [-sectionHalfWidth, 0] },   // left  
      { normal: [0, 1], point: [0, sectionHalfHeight] },    // top
      { normal: [0, -1], point: [0, -sectionHalfHeight] }   // bottom
    ];

    boundaries.forEach(boundary => {
      const [nx, ny] = boundary.normal;
      const [px, py] = boundary.point;

      // Line equation: (x, y) = (offsetX, offsetY) + t * (perpCos, perpSin)
      // Plane equation: nx * (x - px) + ny * (y - py) = 0
      // Solve for t: nx * (offsetX + t * perpCos - px) + ny * (offsetY + t * perpSin - py) = 0

      const denominator = nx * perpCos + ny * perpSin;
      if (Math.abs(denominator) > 0.0001) {
        const t = (nx * (px - offsetX) + ny * (py - offsetY)) / denominator;
        const intersectionX = offsetX + t * perpCos;
        const intersectionY = offsetY + t * perpSin;

        // Check if intersection is within section bounds
        if (intersectionX >= -sectionHalfWidth && intersectionX <= sectionHalfWidth &&
          intersectionY >= -sectionHalfHeight && intersectionY <= sectionHalfHeight) {
          intersections.push({ t, x: intersectionX, y: intersectionY });
        }
      }
    });

    // Sort intersections by t value and take the two extremes
    intersections.sort((a, b) => a.t - b.t);

    if (intersections.length >= 2) {
      const start = [intersections[0].x, intersections[0].y, 0.01];
      const end = [intersections[intersections.length - 1].x, intersections[intersections.length - 1].y, 0.01];
      return [start, end];
    }

    // Fallback to c-based line if no intersections found
    const start = [
      offsetX + perpCos * (-c),
      offsetY + perpSin * (-c),
      0.01
    ];
    const end = [
      offsetX + perpCos * c,
      offsetY + perpSin * c,
      0.01
    ];

    return [start, end];
  }, [naRotation, naPosition, c]);

  // Calculate label positions at the actual extremes of the strain diagram line
  const compressionLabelPos = useMemo(() => {
    if (strainDiagramLine.length >= 2) {
      const [start, end] = strainDiagramLine;

      // Determine which end has more compression (further from neutral axis in negative direction)
      const cosAngle = Math.cos(naRotation);
      const sinAngle = Math.sin(naRotation);
      const perpCos = -sinAngle;
      const perpSin = cosAngle;

      const offsetX = -sinAngle * naPosition;
      const offsetY = cosAngle * naPosition;

      // Calculate distances from neutral axis for both ends
      const startDist = (start[0] - offsetX) * perpCos + (start[1] - offsetY) * perpSin;
      const endDist = (end[0] - offsetX) * perpCos + (end[1] - offsetY) * perpSin;

      // Use the end with more negative distance (compression side)
      const compressionEnd = startDist < endDist ? start : end;

      return [
        compressionEnd[0] - perpCos * 0.3,
        compressionEnd[1] - perpSin * 0.3,
        0.02
      ];
    }
    return [0, 0, 0];
  }, [strainDiagramLine, naRotation, naPosition]);

  const tensionLabelPos = useMemo(() => {
    if (strainDiagramLine.length >= 2) {
      const [start, end] = strainDiagramLine;

      // Determine which end has more tension (further from neutral axis in positive direction)
      const cosAngle = Math.cos(naRotation);
      const sinAngle = Math.sin(naRotation);
      const perpCos = -sinAngle;
      const perpSin = cosAngle;

      const offsetX = -sinAngle * naPosition;
      const offsetY = cosAngle * naPosition;

      // Calculate distances from neutral axis for both ends
      const startDist = (start[0] - offsetX) * perpCos + (start[1] - offsetY) * perpSin;
      const endDist = (end[0] - offsetX) * perpCos + (end[1] - offsetY) * perpSin;

      // Use the end with more positive distance (tension side)
      const tensionEnd = startDist > endDist ? start : end;

      return [
        tensionEnd[0] + perpCos * 0.3,
        tensionEnd[1] + perpSin * 0.3,
        0.02
      ];
    }
    return [0, 0, 0];
  }, [strainDiagramLine, naRotation, naPosition]);

  // Calculate actual strain values at the extreme positions
  const extremeStrains = useMemo(() => {
    if (strainDiagramLine.length >= 2) {
      const [start, end] = strainDiagramLine;

      const cosAngle = Math.cos(naRotation);
      const sinAngle = Math.sin(naRotation);
      const perpCos = -sinAngle;
      const perpSin = cosAngle;

      const offsetX = -sinAngle * naPosition;
      const offsetY = cosAngle * naPosition;

      // Calculate actual distances from neutral axis
      const startDist = (start[0] - offsetX) * perpCos + (start[1] - offsetY) * perpSin;
      const endDist = (end[0] - offsetX) * perpCos + (end[1] - offsetY) * perpSin;

      // Calculate strains: ε = (distance/c) * εmax
      const startStrain = (startDist / c) * maxStrain;
      const endStrain = (endDist / c) * maxStrain;

      const compressionStrain = startDist < endDist ? startStrain : endStrain;
      const tensionStrain = startDist > endDist ? startStrain : endStrain;

      return { compressionStrain, tensionStrain };
    }
    return { compressionStrain: -maxStrain, tensionStrain: maxStrain };
  }, [strainDiagramLine, naRotation, naPosition, c, maxStrain]);

  return (
    <>
      <mesh ref={meshRef} material={material}>
        <planeGeometry args={[4, 4, 64, 64]} />
      </mesh>

      {showNeutralAxis && (<>
        <Line points={neutralAxisLine} color="black" lineWidth={2} dashed dashSize={0.1} gapSize={0.05} />
        <Line points={[compression.x, compression.y, 0.01, intersection.x, intersection.y, 0.01]} color="black" lineWidth={1} />
        <Text position={[(compression.x + intersection.x) / 2, (compression.y + intersection.y) / 2, 0.01]} fontSize={0.15} color="black" anchorX="center" anchorY="middle"          >
          c = {(-c).toFixed(1)}
        </Text>
      </>)}

      {showStrainDiagram && (
        <>
          <Line points={strainDiagramLine} color="white" lineWidth={2} />
          <Text position={compressionLabelPos} fontSize={0.15} color="black" anchorX="center" anchorY="middle"          >
            εc = {(extremeStrains.compressionStrain).toFixed(5)}
          </Text>
          <Text position={tensionLabelPos} fontSize={0.15} color="black" anchorX="center" anchorY="middle"          >
            εt = {(extremeStrains.tensionStrain).toFixed(5)}
          </Text>
        </>
      )}
    </>
  );
}

function Scene() {
  return (
    <Canvas camera={{ position: [0, 0, 5], fov: 50, zoom: 0.5 }} style={{ height: '90vh' }}>
      <ambientLight intensity={0.8} />
      <Sky sunPosition={[0, 1, 0]} />
      <OrbitControls enableDamping dampingFactor={0.05} />
      <NeutralAxisVisualization />
    </Canvas>
  );
}

export default function App() {
  return (
    <div className="w-full h-screen bg-gray-900 relative">
      <Scene />
      <div className="absolute top-4 left-4 text-white">
        <h1 className="text-xl font-bold mb-2">Neutral Axis Strain Visualization</h1>
        <div className="text-sm space-y-1">
          <div>🔵 Blue = Maximum Compression</div>
          <div>🟢 Green = Zero Strain (Neutral)</div>
          <div>🔴 Red = Maximum Tension</div>
          <div>Black Dashed Line = Neutral Axis</div>
          <div>Black Solid Line = c</div>

          <div>White Line = Mid Section Strain</div>
        </div>
      </div>
    </div>
  );
}

/**
 * Calculates the true extreme compression and tension fiber corners
 * for a square centered at (0,0), using the correct perpendicular normal.
 *
 * @param {number} theta     Rotation of the neutral axis (radians CCW from +X).
 * @param {number} d         Offset of the neutral axis along its normal.
 * @param {number} halfSize  Half‐width of your square (e.g. 2 for a 4×4 plane).
 * @returns {{
 *   compression: { x: number, y: number, dist: number },
 *   tension:     { x: number, y: number, dist: number }
 * }}
 */
function getExtremeFibers(theta, d, halfSize) {
  // perpendicular (normal) to the neutral‐axis direction
  const nPerp = {
    x: -Math.sin(theta),
    y: Math.cos(theta)
  };

  // the four corners of the square
  const corners = [
    { x: halfSize, y: halfSize },
    { x: halfSize, y: -halfSize },
    { x: -halfSize, y: halfSize },
    { x: -halfSize, y: -halfSize },
  ];

  // initialize trackers
  let min = { dist: Infinity, x: 0, y: 0 };
  let max = { dist: -Infinity, x: 0, y: 0 };

  // scan each corner
  for (let { x, y } of corners) {
    // signed perpendicular distance from corner to NA
    const dist = x * nPerp.x + y * nPerp.y - d;

    if (dist < min.dist) {
      min = { dist, x, y };
    }
    if (dist > max.dist) {
      max = { dist, x, y };
    }
  }

  return {
    compression: { x: min.x, y: min.y, dist: min.dist },
    tension: { x: max.x, y: max.y, dist: max.dist }
  };
}

/**
 * Calculates the signed perpendicular distance from a point to the neutral axis
 * and returns the intersection point of the perpendicular through that point.
 *
 * Neutral axis is defined by: x·cos(θ) + y·sin(θ) = d
 * We use the same normal as in getExtremeFibers: nPerp = [−sinθ, cosθ].
 *
 * @param {number} x      – x‐coordinate of the point
 * @param {number} y      – y‐coordinate of the point
 * @param {number} theta  – neutral axis rotation in radians (CCW from +X)
 * @param {number} d      – offset of the neutral axis along its normal
 * @returns {{
 *   distance: number,       // signed perpendicular distance
 *   intersection: { x: number, y: number }
 * }}
 */
function getDistanceAndClosestPoint(x, y, theta, d) {
  // normal to the NA direction
  const nPerpX = -Math.sin(theta);
  const nPerpY = Math.cos(theta);

  // signed distance from (x,y) to NA
  const distance = x * nPerpX + y * nPerpY - d;

  // drop a perpendicular from (x,y) to the NA:
  // intersection = (x,y) − distance * nPerp
  const intersectionX = x - distance * nPerpX;
  const intersectionY = y - distance * nPerpY;

  return {
    distance,
    intersection: { x: intersectionX, y: intersectionY }
  };
}

Comments (0)

Loading comments...