Finite Element Method analysis

Question

Please take a look at this class and point out its flaws. I am developing a simple Finite Element Method analysis with PHP.

What this class basically does:

Given an input set (coordinates, properties, plane type), calculate an output/analysis set (area, matrix, etc.).

Possible hints on flaws:

1. This is a full static class
2. Takes array as input parameter, outputs an array too
3. Perhaps the name is too long

If this code is considered a "good design pattern", I will use this similar approach to all other codes in my project.

<?php

class CSTElementStiffnessCalculator {

    public static $sampleInput_StressIsotropic = array(
        'planeType' => 'stress',
        'firstNode' => array(
            'x' => 100,
            'y' => -6.25,
            'id' => 2
        ),
        'secondNode' => array(
            'x' => 200,
            'y' => -12.5,
            'id' => 3
        ),
        'thirdNode' => array(
            'x' => 200,
            'y' => 12.5,
            'id' => 4
        ),
        'property' => array(
            'type' => 'isotropic',
            'E' => 200,
            'V' => 0.3,
            'thickness' => 10
        )
    );
    public static $sampleInput_StrainOrthotropic = array(
        'planeType' => 'strain',
        'firstNode' => array(
            'x' => 0,
            'y' => -6.25,
            'id' => 1
        ),
        'secondNode' => array(
            'x' => 100,
            'y' => -6.25,
            'id' => 5
        ),
        'thirdNode' => array(
            'x' => 100,
            'y' => 6.25,
            'id' => 6
        ),
        'property' => array(
            'type' => 'orthotropic',
            'Ex' => 200,
            'Ey' => 100,
            'Ez' => 50,
            'Vxy' => 0.35,
            'Vxz' => 0.3,
            'Vyz' => 0.25,
            'Gxy' => 50,
            'thickness' => 1
        )
    );

    public static function calculate($inputParams)
    {
        $x1 = $inputParams['firstNode']['x'];
        $y1 = $inputParams['firstNode']['y'];
        $id1 = $inputParams['firstNode']['id'];

        $x2 = $inputParams['secondNode']['x'];
        $y2 = $inputParams['secondNode']['y'];
        $id2 = $inputParams['secondNode']['id'];

        $x3 = $inputParams['thirdNode']['x'];
        $y3 = $inputParams['thirdNode']['y'];
        $id3 = $inputParams['thirdNode']['id'];

        $planeType = $inputParams['planeType'];
        $property = $inputParams['property'];
        $thickness = $property['thickness'];

        /////////////////////////////////////////////////////////////

        $DoubleAMatrix = self::makeDoubleAMatrix(
            $x1, $y1, $x2, $y2, $x3, $y3
        );

        $doubleA = $DoubleAMatrix->det();
        $Area = $doubleA / 2;

        $BMatrix = self::makeBMatrix(
            $Area, $x2-$x1, $x3-$x2, $x1-$x3, $y2-$y1, $y3-$y2, $y1-$y3
        );

        $DMatrix = self::makeDMatrix(
            $planeType, $property
        );

        /////////////////////////////////////////////////////////////

        $BTrans = $BMatrix->transpose();

        $tempMat = $BTrans->times($DMatrix); // [B]T x [D]
        $tempMat = $tempMat->times($BMatrix); // [B]T x [D] x [B]
        $tA = $thickness * $Area;
        $KMatrix = $tempMat->times($tA); // t x A x [B]T x [D] x B

        /////////////////////////////////////////////////////////////

        $stiffnessIds = array();

        $stiffnessIds[] = $id1 * 2 - 1;
        $stiffnessIds[] = $id1 * 2;

        $stiffnessIds[] = $id2 * 2 - 1;
        $stiffnessIds[] = $id2 * 2;

        $stiffnessIds[] = $id3 * 2 - 1;
        $stiffnessIds[] = $id3 * 2;

        /////////////////////////////////////////////////////////////

        return array(
            'd_matrix' => $DMatrix->getArrayCopy(),
            'b_matrix' => $BMatrix->getArrayCopy(),
            'area' => $Area,
            'local_stiff_ids' => $stiffnessIds,
            'local_stiff_matrix' => $KMatrix->getArrayCopy()
        );

    }

    private static function makeDoubleAMatrix($x1, $y1, $x2, $y2, $x3, $y3)
    {
        $data = array(
            array(1, $x1, $y1),
            array(1, $x2, $y2),
            array(1, $x3, $y3),
        );
        return new \Matrix($data);
    }

    private static function makeBMatrix($A, $x12, $x23, $x31, $y12, $y23, $y31)
    {
        $multiplier = 1 / (2 * $A);
        $data = array(
            array( -$y23,     0, -$y31,     0, -$y12,     0 ),
            array(     0,  $x23,     0,  $x31,     0,  $x12 ),
            array(  $x23, -$y23,  $x31, -$y31,  $x12, -$y12 ),
        );
        $matrix = new \Matrix($data);
        return $matrix->times($multiplier);
    }

    private static function makeDMatrix($planeType, $property)
    {
        if ($planeType == 'stress' && $property['type'] == 'isotropic')
        {
            return self::makeDMatrix_Stress_Isotropic(
                $property['E'],
                $property['V']
            );
        }
        if ($planeType == 'stress' && $property['type'] == 'orthotropic')
        {
            return self::makeDMatrix_Stress_Orthotropic(
                $property['Ex'],
                $property['Ey'],
                $property['Vxy'],
                $property['Gxy']
            );
        }
        if ($planeType == 'strain' && $property['type'] == 'isotropic')
        {
            return self::makeDMatrix_Strain_Isotropic(
                $property['E'],
                $property['V']
            );
        }
        if ($planeType == 'strain' && $property['type'] == 'orthotropic')
        {
            return self::makeDMatrix_Strain_Orthotropic(
                $property['Ex'],
                $property['Ey'],
                $property['Ez'],
                $property['Vxy'],
                $property['Vxz'],
                $property['Vyz'],
                $property['Gxy']
            );
        }
    }

    private static function makeDMatrix_Stress_Isotropic($E, $V)
    {
        $multiplier = $E / ( 1 - $V * $V );
        $data = array(
            array(  1, $V, 0 ),
            array( $V,  1, 0 ),
            array(  0,  0, (1-$V)/2 ),
        );
        $matrix = new \Matrix($data);
        return $matrix->times($multiplier);
    }

    ///////////////////////////////////////////////////////////////

    private static function makeDMatrix_Stress_Orthotropic($Ex, $Ey, $Vxy, $Gxy)
    {
        $n = $Ex / $Ey;
        $temp = 1 - $n * $Vxy * $Vxy;
        $v00 = $Ex/$temp;
        $v01 = $v10 = $Vxy*$Ex/$temp;
        $v11 = $Ey/$temp;
        $data = array(
            array( $v00, $v01, 0 ),
            array( $v10, $v11, 0 ),
            array(    0,    0, $Gxy ),
        );
        return new \Matrix($data);
    }

    ///////////////////////////////////////////////////////////////

    private static function makeDMatrix_Strain_Isotropic($E, $V)
    {
        $multiplier = $E / ( ( 1 + $V ) * ( 1 - 2*$V ) );
        $v22 = ( 1 - 2*$V ) / 2;
        $data = array(
            array( 1-$V,   $V, 0 ),
            array(   $V, 1-$V, 0 ),
            array(    0,    0, $v22 ),
        );
        $matrix = new \Matrix($data);
        return $matrix->times($multiplier);
    }

    ///////////////////////////////////////////////////////////////

    private static function makeDMatrix_Strain_Orthotropic($Ex, $Ey, $Ez, $Vxy, $Vxz, $Vyz, $Gxy)
    {
        $D = ( 1 - $Vxz * $Vxz * ($Ex/$Ez) ) * ( 1 - $Vyz * $Vyz * ($Ey/$Ez) ) -
            ( $Vxy + $Vxz * $Vyz * ($Ey/$Ez) ) * ( $Vxy * ($Ex/$Ey) + $Vyz * $Vxz * ($Ex/$Ez));

        $n = 1 / $D;

        $v00 = $n * ( 1 - $Vyz * $Vyz * ($Ey/$Ez) ) * $Ex;
        $v01 = $n * ( $Vxy + $Vxz * $Vyz * ($Ey/$Ez) ) * $Ex;
        $v10 = $n * ( $Vxy * ($Ex/$Ey) + $Vyz * $Vxz * ($Ex/$Ez) ) * $Ey;
        $v11 = $n * ( 1 - $Vxz * $Vxz * ($Ex/$Ez) ) * $Ey;

        $data = array(
            array( $v00, $v01, 0 ),
            array( $v10, $v11, 0 ),
            array(    0,    0, $Gxy ),
        );
        return new \Matrix($data);
    }
}

Answer 1

No, this code is not considered a good design pattern.

It isn't even considered a design pattern. The only thing you did is copy paste a bunch of functions into a class statement and made them static.

Now, since I mostly work on applications not dealing with 'ginite element methods' not everything written down here wil be correct ;) So add a pinch of salt to it!

What do we want to do?

Given 3 Nodes, a PlanetType and Property we need to be able to get the D-matrix B-matrix, area, local stiffness id's and the local stiffness matrix.

So, what do we need?

From the 'what do we want to do' we can extract some objects/entities that we need:

• Node
• PlanetType
• Property
• DMatrix
• BMatric
• Area
• LocalStiffnessIds
• LocalStiffnessMatrix

For each item, we will have to decide what the best way of passing it arround will be. Our node looks like something rather important that is used multiple times in our application. So let's for the sake of it, write an interface for our Node:

<?php namespace FEM\Interfaces;

interface Node {
    /**
     * @param int $id int
     * @param float $x  x
     * @param float $y  y
     */
    public function __construct($id, $x, $y);

    /**
     * @return int id
     */
    public function getId();

    /**
     * @return float x
     */
    public function getX();

    /**
     * @return float y
     */
    public function getY();
}

Now we have our interface ready, let's create a class that implements this interface, so we can atleast start creating our nodes:

<?php namespace FEM;

use FEM\Interfaces\Node as NodeInterface;

class Node implements NodeInterface
{

    /**
     * @var int
     */
    private $id;

    /**
     * @var float
     */
    private $x;

    /**
     * @var float
     */
    private $y;

    public function __construct($id, $x, $y)
    {
        $this->id = (int)$id;   //I'm doing a simple type casting to be sure it's the correct type
        $this->x = $x;          //However, it would be more wise to perform some more explicit checks
        $this->y = $y;          //e.g. $x and $y both 0 wouldn't make a lot of sense I guess?
    }

    public function getId()
    {
        return $this->id;
    }

    public function getX()
    {
        return $this->x;
    }

    public function getY()
    {
        return $this->y;
    }
}

Now we can start creating our nodes:

<?php

use FEM\Node;

$node1 = new Node(2, 100, -6.25);
$node2 = new Node(3, 200, -12.5);
$node3 = new Node(4, 200, -12.5);

$node1->getId(); //2
$node2->getId(); //3
$node3->getId(); //4

We are on fire!

Now that our application knows what Nodes are and can do (interface) we can do the same for our property. In our code, we have 2 different kind of Properties. Isotropic and Orthotropic.

this gives us 2 classes that each extend an AbstractProperty class:

<?php namespace FEM;

use FEM\Interfaces\Property;
use FEM\Interfaces\IsotropicProperty as IsotropicPropertyInterface;
use FEM\Interfaces\OrthotropicProperty as OrthotropicPropertyInterface;

abstract class AbstractProperty implements Property
{

}

class IsotropicProperty extends AbstractProperty implements IsotropicPropertyInterface
{

}

class OrthotropicProperty extends AbstractProperty implements OrthotropicPropertyInterface
{

}

So, know we populate our class with everything itneeds, put all the methods and stuff that all properties have (thickness, ...) in the AbstractProperty. We define our interfaces so that other programmers (read, me in a few weeks) know what every property can do.

Now, the matrix fun. The factory pattern here would be a good approach. You will have a nodeToMatrixFactory (or something like that) with a couple of statis methods:

<?php namespace FEM;

use FEM\Interfaces\Node as NodeInterface;

class NodeToMatrixFactory
{
    public static function createDoubleAMatrix(NodeInterface $node1,NodeInterface $node2,NodeInterface $node3)
    {
        return new \Matrix(array(
            array(1, $node1->getX(), $node1->getY()),
            array(1, $node2->getX(), $node2->getY()),
            array(1, $node3->getX(), $node3->getY())
        ));
    }
}

Note how we can have some nice typehinting there? We tell PHP/ourself that we need 3 objects that implement the Node interface. this means that inside the function we now are 100% sure that $node1->getX() and getY() return a float.

Calculate that piece of $@{!§Node&@('Propertyé"@Stuff

Now our calculate class becomes more of a controller that puts everything together, it creates our Nodes (from a database, file, ...) our properties and ... It then uses the MatrixFactory to get the nice Matrixes for us and then passes it to some other class we wrote to do some more magic, ...

Try to encapsulate every little piece of software in its own little space.

If you can't explain what a class represents/does without using the word 'AND', Chop it into smaller pieces

Answer 2

Someone else will most likely point out several more points to your code.

However, here are two things that have jumped out:

1. Your static class is not scalable. This is one flawed aspect of the class that you have written. Also, OOP means Object Oriented, right? You can't initialize an object out of a static class which in turn makes the class untestable under PHPunit. It also breaks several SOLID design principles. (you can run searches on the site; this has been discussed several times)

2. Some of your functions have too many variables.

   2 ways to shorten this:

   1. use an array or
   2. make an object that holds these variables