main.cpp 2.41 KB
Newer Older
Alexandre Halbach's avatar
Alexandre Halbach committed
1
// This code gives the mechanical vibration eigenfrequencies and eigenmodes of a 
Alexandre Halbach's avatar
Alexandre Halbach committed
2
// 3D disk that is clamped at its outer face.
Alexandre Halbach's avatar
Alexandre Halbach committed
3 4


5
#include "sparselizardbase.h"
6 7 8 9 10 11


using namespace mathop;

void sparselizard(void)
{	
Alexandre Halbach's avatar
Alexandre Halbach committed
12
    // The domain regions as defined in 'disk.geo':
13 14 15
    int vol = 1, sur = 2, top = 3;
    
    // The mesh can be curved!
Alexandre Halbach's avatar
Alexandre Halbach committed
16
    mesh mymesh("disk.msh");
17 18 19 20 21 22 23 24 25 26 27 28 29
    
    // Nodal shape functions 'h1' with 3 components.
    // Field u is the membrane deflection.
    field u("h1xyz");

    // Use interpolation order 3 on 'vol', the whole domain:
    u.setorder(vol, 3);
    
    // Clamp on surface 'sur' (i.e. 0 valued-Dirichlet conditions):
    u.setconstraint(sur);
  
    // E is Young's modulus. nu is Poisson's ratio. rho is the volumic mass.
    parameter E, nu, rho;
Alexandre Halbach's avatar
Alexandre Halbach committed
30
    E|vol = 150e9; nu|vol = 0.3; rho|vol = 2330;
31 32 33 34
  
    formulation elasticity;

    // The linear elasticity formulation is classical and thus predefined:
35
    elasticity += integral(vol, predefinedelasticity(dof(u), tf(u), E, nu));
36
    // Add the inertia terms:
Alexandre Halbach's avatar
Alexandre Halbach committed
37
    elasticity += integral(vol, -rho*dtdt(dof(u))*tf(u));
38 39 40 41 42 43 44

    elasticity.generate();
    
    // Get the stiffness and mass matrix:
    mat K = elasticity.K();
    mat M = elasticity.M();

Alexandre Halbach's avatar
Alexandre Halbach committed
45
    // Remove the rows and columns corresponding to the 0 constraints:
46 47 48
    K.removeconstraints();
    M.removeconstraints();
    
Alexandre Halbach's avatar
Alexandre Halbach committed
49
    // Create the object to solve the generalised eigenvalue problem K*x = lambda*M*x :
50 51 52 53 54
    eigenvalue eig(K, M);
    
    // Compute the 10 eigenvalues closest to the target magnitude 0.0 (i.e. the 10 first ones):
    eig.compute(10, 0.0);
    
Alexandre Halbach's avatar
Alexandre Halbach committed
55 56
    // Print the eigenfrequencies:
    eig.printeigenfrequencies();
57 58 59 60 61 62 63 64
    
    // The eigenvectors are real thus we only need the real part:
    std::vector<vec> myeigenvectors = eig.geteigenvectorrealpart();
    
    // Loop on all eigenvectors found:
    for (int i = 0; i < myeigenvectors.size(); i++)
    {
    	// Transfer the data from the ith eigenvector to field u:
Alexandre Halbach's avatar
Alexandre Halbach committed
65
        u.setdata(top, myeigenvectors[i]);
Alexandre Halbach's avatar
Alexandre Halbach committed
66
        // Write the deflection on the top surface of the membrane with an order 3 interpolation:
67 68
        u.write(top, "u"+std::to_string(i)+".pos", 3);
    }
69 70
    
    // Code validation line. Can be removed.
Alexandre Halbach's avatar
Alexandre Halbach committed
71
    std::cout << (eig.geteigenvaluerealpart()[0] < 6.25240e+06 && eig.geteigenvaluerealpart()[0] > 6.25235e+06);
72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92
}

int main(void)
{	
    SlepcInitialize(0,{},0,0);

    sparselizard();

    SlepcFinalize();

    return 0;
}