Haplo

Le but de cet exemple est de valider les concepts mis en oeuvre dans Haplo. Evidement, le problème choisi est extrêmement simple. Cependant il met en oeuvre les principes suivants:

double dualisation des conditions aux limites,
numérotation des inconnues (déplacements et coefficients de lagrange),
matrice en stockage ligne de ciel,
inversion d'un système linéaire avec un algorithme de type LDLT optimisé,
et bien-sûr les mécanismes de modules et d'objets d'Haplo.

Description du problème

La figure suivante illustre le problème d'élasticité à résoudre.

On cherche le déplacement du noeud sur lequel la force est appliquée. Ce résultat est DX=8

Utilisation d'Haplo pour le résoudre

< H A P L O > Version 0.7.1 (En) Copyright (C) 1998-2003 Nicolas LAURENT	message de bienvenue.
>> use("fem"); use: <F.E.M for Mechanics> version 0.7.1 (c) Nicolas LAURENT fem: Définition de l'élément 'ME1DL2'	On charge le module de calcul par éléments finis.
>> ma=meshing(1); Maillage 1D (0 noeuds, 0 mailles) >> N=4; 4 >> node_add(ma, [0]); >> for(i = 1; i <= N; i = i+1) >> { >> node_add(ma, [2*i]); >> mesh_add(ma, [i, i+1]); >> } >> print(ma); +----------+ \| MAILLAGE \| +----------+ 5 nodes ------- 1. 0.000 2. 2.000 3. 4.000 4. 6.000 5. 8.000 4 meshes -------- 1. 1 2 2. 2 3 3. 3 4 4. 4 5	Création du maillage et affichage.
>> mo=model(ma); Model (Mecanic 1D) >> print(mo); +-------+ \| MODEL \| +-------+ 4 elements ---------- 1. ME1DL2: N1 N2 2. ME1DL2: N2 N3 3. ME1DL2: N3 N4 4. ME1DL2: N4 N5	Création d'un modèle mathématique s'appuyant sur ce maillage. Le choix des élements n'est pas spécifié dans cet exemple: le moteur de calcul choisira des élements finis compatibles avec les mailles.
>> mat=material(); material >> set(mat, "E", 1); >> print(mat); +----------+ \| Material \| +----------+ E = 1 >> set(mo, mat);	Définition d'un matériau et affectation sur tous les élements du modèle.
>> pro=property(); property >> set(pro, "S", 1); >> print(pro); +----------+ \| Property \| +----------+ S = 1 >> set(mo, pro);	Définition d'une propriété géométrique et affectation sur tous les élements du modèle.
>> bc=conditions(); Boundry conditions set (0 constraint) >> set(bc, 1, "DX", 0); >> print(bc); +---------------------+ \| Boundary conditions \| +---------------------+ 1 condition ----------- DX1 = 0	Définition de conditions aux limites. Pour l'instant, celles-ci ne sont pas affectées au modèle.
>> load=loading(); Loading (0 loads) >> add(load, N+1, "FX", 1.0); >> print(load); +---------+ \| Loading \| +---------+ 1 load ------ FX5 = 1.000000e+00	Définition d'un chargement. Pour l'instant, celui-ci n'est pas affecté au modèle.
>> P=problem(mo, bc, load); problem: Numbering: <0.01s Problem (5 dof, 1 boundary condition)	Création d'un problème mécanique complet: modèle mathématique+conditions aux limites+ chargement. A partir ce cette instruction, les inconnues sont numérotées.
>> K=stiffness(P); stiffness: BC coefficient: alpha=7.500000e-01 stiffness: Assembly: <0.01s Matrix (7x7, 53.6%) >> profile(K); ( X ) ( XX ) ( XXX ) ( XXX ) ( XX ) ( XX ) ( XX ) >> print(K); ( -7.50e-01 ) ( +7.50e-01 +5.00e-01 ) ( +7.50e-01 +7.50e-01 -7.50e-01 ) ( -5.00e-01 +0.00e+00 +1.00e+00 ) ( -5.00e-01 +1.00e+00 ) ( -5.00e-01 +1.00e+00 ) ( -5.00e-01 +5.00e-01 )	Calcul de la matrice de rigidité du problème. Les inconnues nommées LAMBDA1 et LAMBDA2 correspondent à la double dualisation des conditions aux limites.
>> for(i=1; i <= dof(P); i=i+1) >> print(i + " <-> " + dof(P, i)); 1 <-> LAMBDA1_1 2 <-> DX1 3 <-> LAMBDA2_1 4 <-> DX2 5 <-> DX3 6 <-> DX4 7 <-> DX5	Affichage des inconnues et de leurs numérotations.
>> F=vector(P); Vector (7) >> print(F); [ +0.00e+00 ] [ +0.00e+00 ] [ +0.00e+00 ] [ +0.00e+00 ] [ +0.00e+00 ] [ +0.00e+00 ] [ +1.00e+00 ]	Calcul du vecteur de chargement (second membre).
>> U=solve_ldlt(K, F); solve_ldlt: Factorization: <0.01s solve_ldlt: Inversion: <0.01s Vector (7)	Dans le cas de l'élasticité linéaire, il s'agit d'inverser le système. Dans cet exemple, on choisit un solveur direct (le solveur itératif est également disponible).
>> dx=get(U, dof(P, N+1, "DX")); 8	On affiche les résultats recherchés. On obtient effectivement un déplacement de 8!
>> quit; CPU total = 0.03s user = 0.03s (79.4%) sys = 0.01s (20.6%) Number of malloc()....: 959 Number of free()......: 959 Number of realloc()...: 0 Maximal RSS...........: 494.41 kb No memory leaks. End of processing.	Messages récapitulatifs sur l'utilisation des ressources systèmes.

Dernière modification: $Date: 2004/06/20 15:04:59 $

This page conforms to xhtml 1.1 and css2 standards