Glassvision/JuliaFEM.jl
1
0
Fork 0
mirror of https://github.com/JuliaFEM/JuliaFEM.jl.git synced 2026-08-06 04:21:33 +00:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'master' of https://github.com/JuliaFEM/JuliaFEM.jl

Browse Source
This commit is contained in:
ovainola
2015-08-25 21:09:11 +03:00
parent f5137a2603 9e7280c53b
commit 0b77403324
18 changed files with 3821 additions and 1960 deletions
Download Patch File Download Diff File Expand all files Collapse all files
.gitignore
+2
View File
@@ -1,3 +1,5 @@
*~ *~
.DS_Store .DS_Store
.ipynb_checkpoints .ipynb_checkpoints
docs/build/html
*.swp
CONTRIBUTING.rst
+29
View File
@@ -9,6 +9,35 @@ For now, read
https://github.com/JuliaLang/julia/blob/master/CONTRIBUTING.md https://github.com/JuliaLang/julia/blob/master/CONTRIBUTING.md
How to contribute
-----------------
Here are the basic steps for contributing to JuliaFEM:
1) Create an account or sign in to `GitHub <https://github.com/>`_.
2) Go to `Git home page <http://git-scm.com/>`_ and download the Git installer. Run the installer to get Git on your computer. It is a version control system used by GitHub. To learn its basics, go through this `Git tutorial <https://try.github.io/levels/1/challenges/1>`_.
3) Install Julia (v0.4+) to your computer. At `Julia readme
<https://github.com/JuliaLang/julia/blob/master/README.md>`_ you'll find complete instructions for installing it for your platform.
4) Go to the `JuliaFEM GitHub page <https://github.com/JuliaFEM/JuliaFEM.jl>`_. At the top-right corner, press the ``Fork``-button to fork your own copy of JuliaFEM to your repository.
5) Clone JuliaFEM from your repository to your computer. Navigate to the folder you want to clone it to, and type the following command (inserting your GitHub username to its place):
``git clone https://github.com/your_github_username/JuliaFEM.jl.git``
6) You can now navigate to JuliaFEM in the folder you chose at step 5. There you'll find the same contents as you see in your GitHub JuliaFEM repository. Now, locate the file you want to modify, open it with your desired text editor, make the changes and save the new version. If you type ``git status``, you'll see that the files you've created or modified are listed under ``untracked files``.
7) Add the files you want to update to the staging area by typing ``git add <file1> <file2>...``. If you type ``git status``, you'll see that the files added to the staging area are listed under ``Changes to be committed``. This process also supports wildcard symbols. If you want to add all the files to the staging area, just type ``git add .``. If you want to remove a file from the staging area, type ``git reset <file>``.
8) To store the staged files, commit the files to your repository and add a description message by typing ``git commit -m "your_message_here"``. The message should describe the changes that were made.
9) When you are happy with the commits and want to update them to your repository, type ``git push origin master``.
10) Go to your GitHub JuliaFEM repository. You'll notice that the commit you have made and pushed is now visible above the JuliaFEM file branch. If you click the ``latest commit`` link, you can see the changes made to the file. Finally, click ``Pull request`` to create a pull request of the commits you've made, so that other contributors can review it.
11) If other contributors ask you to make changes to your pull request, just repeat steps 6-9. Your commits will be updated to your original pull request. Do this until everyone is satisfied and your pull request can be merged to the master branch.
Developing Developing
---------- ----------
```bash ```bash
docs/build_notebooks.jl
+1 -1
View File
@@ -55,7 +55,7 @@ function run_notebooks()
# port = 34211+k # we're having some weird port issue with zmq # port = 34211+k # we're having some weird port issue with zmq
# k += 1 # k += 1
try try
run(`runipy -o tutorials/$ipynb --kernel=julia-0.4`) run(`timeout 180 runipy -o tutorials/$ipynb --kernel=julia-0.4`)
status = 0 status = 0
catch catch
println("did not work") println("did not work")
docs/links.md
+3 -1
View File
@@ -22,4 +22,6 @@ Integration
- http://arxiv.org/pdf/1411.1341.pdf - http://arxiv.org/pdf/1411.1341.pdf
Hierarchial shape functions Hierarchial shape functions
- https://www.math.vt.edu/people/adjerids/research/papers/basis.pdf - https://www.math.vt.edu/people/adjerids/research/papers/basis.pdf
- edited by Ari
notebooks/2015-06-25-elasticity-solver-example.ipynb
+1611 -1733
View File
File diff suppressed because one or more lines are too long
notebooks/2015-06-25-shape-functions.ipynb
+440 -46
View File
@@ -1,5 +1,16 @@
{ {
"cells": [ "cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Shape function and integration points\n",
"\n",
"Author(s): Jukka Aho\n",
"\n",
"**Abstract**: Shape functions and element descriptions used in JuliaFEM."
]
},
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 1, "execution_count": 1,
@@ -8,18 +19,325 @@
}, },
"outputs": [], "outputs": [],
"source": [ "source": [
"from sympy import *" "from sympy import *\n",
"#init_printing()"
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 2, "execution_count": 2,
"metadata": { "metadata": {
"collapsed": true "collapsed": false
}, },
"outputs": [], "outputs": [],
"source": [ "source": [
"xi = DeferredVector(\"xi\")" "xi = DeferredVector(r\"xi\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 1D shape function"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Linear 2-node segment (Lagrange family)\n",
"\n",
"| | $\\xi_1$ |\n",
"| ----- | -------:|\n",
"| $N_1$ | -1 |\n",
"| $N_2$ | 1 |\n",
"\n",
"\\begin{equation}\n",
" \\left(\\mathbf{P}\\boldsymbol{\\alpha}\\right)\\left(\\xi_1\\right) = \\alpha_{1}+\\alpha_{2}\\xi_{1}\n",
"\\end{equation}"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"(Matrix([\n",
" [-xi[1]/2 + 1/2],\n",
" [ xi[1]/2 + 1/2]]), Matrix([\n",
" [-1/2],\n",
" [ 1/2]]))"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = Matrix([[1, -1], [1, 1]])\n",
"P = Matrix([1, xi[1]]).T\n",
"N = (P*A.inv()).T\n",
"dN = Matrix([N.diff(xi[1]).T]).T\n",
"N, dN"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Quadratic 3-node segment (Lagrange family)\n",
"\n",
"| | $\\xi_1$ |\n",
"| ----- | -------:|\n",
"| $N_1$ | -1 |\n",
"| $N_2$ | 1 |\n",
"| $N_3$ | 0 |\n",
"\n",
"\\begin{equation}\n",
" \\left(\\mathbf{P}\\boldsymbol{\\alpha}\\right)\\left(\\xi_1\\right) = \\alpha_1 + \\alpha_2\\xi_1 + \\alpha_3\\xi_1^2\n",
"\\end{equation}"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"(Matrix([\n",
" [xi[1]**2/2 - xi[1]/2],\n",
" [xi[1]**2/2 + xi[1]/2],\n",
" [ -xi[1]**2 + 1]]), Matrix([\n",
" [xi[1] - 1/2],\n",
" [xi[1] + 1/2],\n",
" [ -2*xi[1]]]))"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = Matrix([[1, -1, (-1)**2],\n",
" [1, 1, 1**2],\n",
" [1, 0, 0**2]])\n",
"P = Matrix([1, xi[1], xi[1]**2]).T\n",
"N = (P*A.inv()).T\n",
"dN = Matrix([N.diff(xi[1]).T]).T\n",
"N, dN"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### P-elements"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 2D shape functions"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Linear triangle\n",
"\n",
"| | $\\xi_1$ | $\\xi_2$ |\n",
"| ----- | -------:| -------:|\n",
"| $N_1$ | 0 | 0 |\n",
"| $N_2$ | 1 | 0 |\n",
"| $N_3$ | 0 | 1 |"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"(Matrix([\n",
" [-xi[1] - xi[2] + 1],\n",
" [ xi[1]],\n",
" [ xi[2]]]), Matrix([\n",
" [-1, -1],\n",
" [ 1, 0],\n",
" [ 0, 1]]))"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = Matrix([[1, 0, 0], [1, 1, 0], [1, 0, 1]])\n",
"P = Matrix([1, xi[1], xi[2]]).T\n",
"N = (P*A.inv()).T\n",
"dN = Matrix([N.diff(xi[1]).T, N.diff(xi[2]).T]).T\n",
"N, dN"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Quadratic triangle\n",
"\n",
"| | $\\xi_1$ | $\\xi_2$ |\n",
"| ----- | -------:| -------:|\n",
"| $N_1$ | 0 | 0 |\n",
"| $N_2$ | 1 | 0 |\n",
"| $N_3$ | 0 | 1 |\n",
"| $N_4$ | 1/2 | 0 |\n",
"| $N_5$ | 1/2 | 1/2 |\n",
"| $N_6$ | 0 | 1/2 |"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"Matrix([\n",
"[1, 0, 0, 0, 0, 0],\n",
"[1, 1, 0, 1, 0, 0],\n",
"[1, 0, 1, 0, 1, 0],\n",
"[1, 1/2, 0, 1/4, 0, 0],\n",
"[1, 1/2, 1/2, 1/4, 1/4, 1/4],\n",
"[1, 0, 1/2, 0, 1/4, 0]])"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"P = Matrix([1, xi[1], xi[2], xi[1]**2, xi[2]**2, xi[1]*xi[2]]).T\n",
"A = Matrix([\n",
" P.subs({xi[1]: 0, xi[2]: 0}),\n",
" P.subs({xi[1]: 1, xi[2]: 0}),\n",
" P.subs({xi[1]: 0, xi[2]: 1}),\n",
" P.subs({xi[1]: Rational(1,2), xi[2]: 0}),\n",
" P.subs({xi[1]: Rational(1,2), xi[2]: Rational(1,2)}),\n",
" P.subs({xi[1]: 0, xi[2]: Rational(1,2)}),\n",
" ])\n",
"A"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"(Matrix([\n",
" [2*xi[1]**2 + 4*xi[1]*xi[2] - 3*xi[1] + 2*xi[2]**2 - 3*xi[2] + 1],\n",
" [ 2*xi[1]**2 - xi[1]],\n",
" [ 2*xi[2]**2 - xi[2]],\n",
" [ -4*xi[1]**2 - 4*xi[1]*xi[2] + 4*xi[1]],\n",
" [ 4*xi[1]*xi[2]],\n",
" [ -4*xi[1]*xi[2] - 4*xi[2]**2 + 4*xi[2]]]), Matrix([\n",
" [ 4*xi[1] + 4*xi[2] - 3, 4*xi[1] + 4*xi[2] - 3],\n",
" [ 4*xi[1] - 1, 0],\n",
" [ 0, 4*xi[2] - 1],\n",
" [-8*xi[1] - 4*xi[2] + 4, -4*xi[1]],\n",
" [ 4*xi[2], 4*xi[1]],\n",
" [ -4*xi[2], -4*xi[1] - 8*xi[2] + 4]]))"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"N = (P*A.inv()).T\n",
"dN = Matrix([N.diff(xi[1]).T, N.diff(xi[2]).T]).T\n",
"N, dN"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 3D shape functions"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Linear tetrahedra, **tet4**\n",
"\n",
"| | $\\xi_1$ | $\\xi_2$ | $\\xi_2$ |\n",
"| ----- | -------:| -------:| -------:|\n",
"| $N_1$ | 0 | 0 | 0 |\n",
"| $N_2$ | 1 | 0 | 0 |\n",
"| $N_3$ | 0 | 1 | 0 |\n",
"| $N_4$ | 0 | 0 | 1 |"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"Matrix([\n",
"[1, 0, 0, 0],\n",
"[1, 1, 0, 0],\n",
"[1, 0, 1, 0],\n",
"[1, 0, 0, 1]])"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"P = Matrix([1, xi[1], xi[2], xi[3]]).T\n",
"A = Matrix([\n",
" P.subs({xi[1]: 0, xi[2]: 0, xi[3]: 0}),\n",
" P.subs({xi[1]: 1, xi[2]: 0, xi[3]: 0}),\n",
" P.subs({xi[1]: 0, xi[2]: 1, xi[3]: 0}),\n",
" P.subs({xi[1]: 0, xi[2]: 0, xi[3]: 1}),\n",
" ])\n",
"A"
] ]
}, },
{ {
@@ -32,17 +350,15 @@
{ {
"data": { "data": {
"text/plain": [ "text/plain": [
"Matrix([\n", "(Matrix([\n",
"[ -xi[1] - xi[2] - xi[3] + 1],\n", " [-xi[1] - xi[2] - xi[3] + 1],\n",
"[ xi[1]],\n", " [ xi[1]],\n",
"[ xi[2]],\n", " [ xi[2]],\n",
"[ xi[3]],\n", " [ xi[3]]]), Matrix([\n",
"[4*xi[1]*(-xi[1] - xi[2] - xi[3] + 1)],\n", " [-1, -1, -1],\n",
"[ 4*xi[1]*xi[2]],\n", " [ 1, 0, 0],\n",
"[4*xi[2]*(-xi[1] - xi[2] - xi[3] + 1)],\n", " [ 0, 1, 0],\n",
"[4*xi[3]*(-xi[1] - xi[2] - xi[3] + 1)],\n", " [ 0, 0, 1]]))"
"[ 4*xi[1]*xi[3]],\n",
"[ 4*xi[2]*xi[3]]])"
] ]
}, },
"execution_count": 9, "execution_count": 9,
@@ -51,23 +367,29 @@
} }
], ],
"source": [ "source": [
"def c3d10():\n", "N = (P*A.inv()).T\n",
" N1 = 1 - xi[1] - xi[2] - xi[3]\n", "dN = Matrix([N.diff(xi[1]).T, N.diff(xi[2]).T, N.diff(xi[3]).T]).T\n",
" N2 = xi[1]\n", "N, dN"
" N3 = xi[2]\n", ]
" N4 = xi[3]\n", },
" N5 = 4*xi[1]*(1-xi[1]-xi[2]-xi[3])\n", {
" N6 = 4*xi[1]*xi[2]\n", "cell_type": "markdown",
" N7 = 4*xi[2]*(1-xi[1]-xi[2]-xi[3])\n", "metadata": {},
" N8 = 4*xi[3]*(1-xi[1]-xi[2]-xi[3])\n", "source": [
" N9 = 4*xi[1]*xi[3]\n", "### Quadratic Lagrange tetrahedral element, 10 nodes, **tet10**\n",
" N10 = 4*xi[2]*xi[3]\n",
" N = Matrix([N1, N2, N3, N4, N5, N6, N7, N8, N9, N10])\n",
" dN = Matrix([N.diff(xi[1]).T, N.diff(xi[2]).T, N.diff(xi[3]).T]).T\n",
" return N, dN\n",
"\n", "\n",
"N, dN = c3d10()\n", "| | $\\xi_1$ | $\\xi_2$ | $\\xi_2$ |\n",
"N" "| ----- | -------:| -------:| -------:|\n",
"| $N_1$ | 0 | 0 | 0 |\n",
"| $N_2$ | 1 | 0 | 0 |\n",
"| $N_3$ | 0 | 1 | 0 |\n",
"| $N_4$ | 0 | 0 | 1 |\n",
"| $N_5$ | 1/2 | 0 | 0 |\n",
"| $N_6$ | 1/2 | 1/2 | 0 |\n",
"| $N_7$ | 0 | 1/2 | 0 |\n",
"| $N_8$ | 0 | 0 | 1/2 |\n",
"| $N_9$ | 1/2 | 0 | 1/2 |\n",
"| $N_{10}$ | 0 | 1/2 | 1/2 |"
] ]
}, },
{ {
@@ -81,16 +403,16 @@
"data": { "data": {
"text/plain": [ "text/plain": [
"Matrix([\n", "Matrix([\n",
"[ -1, -1, -1],\n", "[1, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n",
"[ 1, 0, 0],\n", "[1, 1, 0, 0, 1, 0, 0, 0, 0, 0],\n",
"[ 0, 1, 0],\n", "[1, 0, 1, 0, 0, 0, 1, 0, 0, 0],\n",
"[ 0, 0, 1],\n", "[1, 0, 0, 1, 0, 0, 0, 0, 1, 0],\n",
"[-8*xi[1] - 4*xi[2] - 4*xi[3] + 4, -4*xi[1], -4*xi[1]],\n", "[1, 1/2, 0, 0, 1/4, 0, 0, 0, 0, 0],\n",
"[ 4*xi[2], 4*xi[1], 0],\n", "[1, 1/2, 1/2, 0, 1/4, 1/4, 1/4, 0, 0, 0],\n",
"[ -4*xi[2], -4*xi[1] - 8*xi[2] - 4*xi[3] + 4, -4*xi[2]],\n", "[1, 0, 1/2, 0, 0, 0, 1/4, 0, 0, 0],\n",
"[ -4*xi[3], -4*xi[3], -4*xi[1] - 4*xi[2] - 8*xi[3] + 4],\n", "[1, 0, 0, 1/2, 0, 0, 0, 0, 1/4, 0],\n",
"[ 4*xi[3], 0, 4*xi[1]],\n", "[1, 1/2, 0, 1/2, 1/4, 0, 0, 0, 1/4, 1/4],\n",
"[ 0, 4*xi[3], 4*xi[2]]])" "[1, 0, 1/2, 1/2, 0, 0, 1/4, 1/4, 1/4, 0]])"
] ]
}, },
"execution_count": 10, "execution_count": 10,
@@ -99,17 +421,89 @@
} }
], ],
"source": [ "source": [
"dN" "P = Matrix([1, xi[1], xi[2], xi[3], xi[1]**2, xi[1]*xi[2], xi[2]**2, xi[2]*xi[3], xi[3]**2, xi[1]*xi[3]]).T\n",
"A = Matrix([\n",
" P.subs({xi[1]: 0, xi[2]: 0, xi[3]: 0}),\n",
" P.subs({xi[1]: 1, xi[2]: 0, xi[3]: 0}),\n",
" P.subs({xi[1]: 0, xi[2]: 1, xi[3]: 0}),\n",
" P.subs({xi[1]: 0, xi[2]: 0, xi[3]: 1}),\n",
"\n",
" P.subs({xi[1]: Rational(1,2), xi[2]: 0, xi[3]: 0}),\n",
" P.subs({xi[1]: Rational(1,2), xi[2]: Rational(1,2), xi[3]: 0}),\n",
" P.subs({xi[1]: 0, xi[2]: Rational(1,2), xi[3]: 0}),\n",
"\n",
" P.subs({xi[1]: 0, xi[2]: 0, xi[3]: Rational(1,2)}),\n",
" P.subs({xi[1]: Rational(1,2), xi[2]: 0, xi[3]: Rational(1,2)}),\n",
" P.subs({xi[1]: 0, xi[2]: Rational(1,2), xi[3]: Rational(1,2)}),\n",
" ])\n",
"A"
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": 11,
"metadata": { "metadata": {
"collapsed": true "collapsed": false
}, },
"outputs": [], "outputs": [
"source": [] {
"data": {
"text/plain": [
"Matrix([\n",
"[(xi[1] + xi[2] + xi[3] - 1)*(2*xi[1] + 2*xi[2] + 2*xi[3] - 1)],\n",
"[ xi[1]*(2*xi[1] - 1)],\n",
"[ xi[2]*(2*xi[2] - 1)],\n",
"[ xi[3]*(2*xi[3] - 1)],\n",
"[ -4*xi[1]*(xi[1] + xi[2] + xi[3] - 1)],\n",
"[ 4*xi[1]*xi[2]],\n",
"[ -4*xi[2]*(xi[1] + xi[2] + xi[3] - 1)],\n",
"[ -4*xi[3]*(xi[1] + xi[2] + xi[3] - 1)],\n",
"[ 4*xi[1]*xi[3]],\n",
"[ 4*xi[2]*xi[3]]])"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"N = (P*A.inv()).T\n",
"dN = Matrix([N.diff(xi[1]).T, N.diff(xi[2]).T, N.diff(xi[3]).T]).T\n",
"factor(N)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"Matrix([\n",
"[ 4*xi[1] + 4*xi[2] + 4*xi[3] - 3, 4*xi[1] + 4*xi[2] + 4*xi[3] - 3, 4*xi[1] + 4*xi[2] + 4*xi[3] - 3],\n",
"[ 4*xi[1] - 1, 0, 0],\n",
"[ 0, 4*xi[2] - 1, 0],\n",
"[ 0, 0, 4*xi[3] - 1],\n",
"[-4*(2*xi[1] + xi[2] + xi[3] - 1), -4*xi[1], -4*xi[1]],\n",
"[ 4*xi[2], 4*xi[1], 0],\n",
"[ -4*xi[2], -4*(xi[1] + 2*xi[2] + xi[3] - 1), -4*xi[2]],\n",
"[ -4*xi[3], -4*xi[3], -4*(xi[1] + xi[2] + 2*xi[3] - 1)],\n",
"[ 4*xi[3], 0, 4*xi[1]],\n",
"[ 0, 4*xi[3], 4*xi[2]]])"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"factor(dN)"
]
} }
], ],
"metadata": { "metadata": {
@@ -128,7 +522,7 @@
"name": "python", "name": "python",
"nbconvert_exporter": "python", "nbconvert_exporter": "python",
"pygments_lexer": "ipython2", "pygments_lexer": "ipython2",
"version": "2.7.9" "version": "2.7.10"
} }
}, },
"nbformat": 4, "nbformat": 4,
notebooks/2015-08-15-interpolation-integration-and-linearization-strategies.ipynb
+548
View File
@@ -0,0 +1,548 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Interpolation and integration algorithms\n",
"\n",
"Author(s): Jukka Aho\n",
"\n",
"**Abstract**: Some strategies to implement automatic differentiation. The number of different choises is caused by a fact that the linearization of function can be done before integration or vice versa, and functions can return values or do in-place modifications. There is probably performance differences between different strategies, but all of them should work."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"using JuliaFEM\n",
"using ForwardDiff"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\"Old good\" elasticity force equilibrium equation $R = T - F$"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"calc_residual_vector_integrand (generic function with 1 method)"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function calc_residual_vector_integrand(el::JuliaFEM.Element, xi)\n",
" # Calculate dN/dX and interpolate material parameters\n",
" dbasisdX = JuliaFEM.get_dbasisdX(el, xi)\n",
" u = el.attributes[\"displacement\"]\n",
" lambda = JuliaFEM.interpolate(el, \"lambda\", xi)\n",
" mu = JuliaFEM.interpolate(el, \"mu\", xi)\n",
"\n",
" # Calculate residual force vector R = T - F\n",
" gradu = u*dbasisdX\n",
" F = I + gradu\n",
" E = 1/2*(gradu' + gradu + gradu'*gradu)\n",
" S = lambda*trace(E)*I + 2*mu*E\n",
" P = F*S\n",
" T = P*dbasisdX'\n",
" return T\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Test case, already well known 2d elasticity in [0,10] x [0,1] grid."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"basis(xi) = [\n",
" (1-xi[1])*(1-xi[2])/4\n",
" (1+xi[1])*(1-xi[2])/4\n",
" (1+xi[1])*(1+xi[2])/4\n",
" (1-xi[1])*(1+xi[2])/4]\n",
"dbasis(xi) = [-(1-xi[2])/4.0 -(1-xi[1])/4.0\n",
" (1-xi[2])/4.0 -(1+xi[1])/4.0\n",
" (1+xi[2])/4.0 (1+xi[1])/4.0\n",
" -(1+xi[2])/4.0 (1-xi[1])/4.0]\n",
"ipoints = 1/sqrt(3)*[-1 -1; 1 -1; 1 1; -1 1]'\n",
"iweights = [1, 1, 1, 1]\n",
"attributes = Dict()\n",
"e = JuliaFEM.Element(1, [1, 2, 3, 4], basis, dbasis, attributes, ipoints, iweights)\n",
"\n",
"E = 90\n",
"nu = 0.25\n",
"mu = E/(2*(1+nu))\n",
"la = E*nu/((1+nu)*(1-2*nu))\n",
"la = 2*la*mu/(la + 2*mu)\n",
"\n",
"e.attributes[\"coordinates\"] = [0.0 0.0; 10.0 0.0; 10.0 1.0; 0.0 1.0]'\n",
"e.attributes[\"lambda\"] = la\n",
"e.attributes[\"mu\"] = mu\n",
"e.attributes[\"displacement\"] = [0.0 0.0; 0.0 0.0; 0.5 0.0; 0.0 0.0]'\n",
"e.attributes[\"displacement nodal force\"] = zeros(2, 4)\n",
"e.attributes[\"displacement tangent stiffness\"] = zeros(8, 8);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Integration\n",
"\n",
"1. take element and function and return value\n",
"2. take function and return function which can be integrated by passing element as a function\n",
"3. do in-place integration, save values to target"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"integrate! (generic function with 1 method)"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function integrate(f::Function, el::JuliaFEM.Element)\n",
" target = []\n",
" for m = 1:length(el.iweights)\n",
" w = el.iweights[m]\n",
" xi = el.ipoints[:, m]\n",
" J = JuliaFEM.interpolate(el, \"coordinates\", xi; derivative=true)\n",
" push!(target, w*f(el, xi)*det(J))\n",
" end\n",
" return sum(target)\n",
"end\n",
"\n",
"function integrate(f::Function)\n",
" function integrate(el::JuliaFEM.Element)\n",
" target = []\n",
" for m = 1:length(el.iweights)\n",
" w = el.iweights[m]\n",
" xi = el.ipoints[:, m]\n",
" J = JuliaFEM.interpolate(el, \"coordinates\", xi; derivative=true)\n",
" push!(target, w*f(el, xi)*det(J))\n",
" end\n",
" return sum(target)\n",
" end\n",
" return integrate\n",
"end\n",
"\n",
"function integrate!(f::Function, el::JuliaFEM.Element, target)\n",
" # set target to zero\n",
" el.attributes[target][:] = 0.0\n",
" for m = 1:length(el.iweights)\n",
" w = el.iweights[m]\n",
" xi = el.ipoints[:, m]\n",
" J = JuliaFEM.interpolate(el, \"coordinates\", xi; derivative=true)\n",
" el.attributes[target][:,:] += w*f(el, xi)*det(J)\n",
" end\n",
"end\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"2x4 Array{Float64,2}:\n",
" -38.2303 -72.8697 79.4912 31.6088\n",
" -17.625 -28.475 37.7 8.4 "
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"integrate(calc_residual_vector_integrand, e)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"2x4 Array{Float64,2}:\n",
" -38.2303 -72.8697 79.4912 31.6088\n",
" -17.625 -28.475 37.7 8.4 "
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"calc_residual_vector = integrate(calc_residual_vector_integrand)\n",
"calc_residual_vector(e)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"2x4 Array{Float64,2}:\n",
" -38.2303 -72.8697 79.4912 31.6088\n",
" -17.625 -28.475 37.7 8.4 "
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"integrate!(calc_residual_vector_integrand, e, \"displacement nodal force\")\n",
"e.attributes[\"displacement nodal force\"]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Linearization\n",
"\n",
"1. take function, element and field, and return partial derivative\n",
"2. take function and field, return function which takes element as argument\n",
"3. do in-place linearization to target, requires function which takes element as argument\n",
"\n",
"In general linearization can be done before integration and vice versa, i.e.\n",
"\n",
" integrate(linearize(f, \"u\"))(e) <-> linearize(integrate(f), \"u\")(e)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"linearize! (generic function with 1 method)"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function linearize(f::Function, el::JuliaFEM.Element, field::ASCIIString)\n",
" dim, nnodes = size(el.attributes[field])\n",
" function helper!(x, y)\n",
" orig = copy(el.attributes[field])\n",
" el.attributes[field] = reshape(x, dim, nnodes)\n",
" y[:] = f(el)\n",
" el.attributes[field] = copy(orig)\n",
" end\n",
" jac = ForwardDiff.forwarddiff_jacobian(helper!, Float64, fadtype=:dual, n=dim*nnodes, m=dim*nnodes)\n",
" return jac(el.attributes[field][:])\n",
"end\n",
"\n",
"function linearize(f::Function, field::ASCIIString)\n",
" function jacobian(el::JuliaFEM.Element, args...)\n",
" dim, nnodes = size(el.attributes[field])\n",
" function helper!(x, y)\n",
" orig = copy(el.attributes[field])\n",
" el.attributes[field] = reshape(x, dim, nnodes)\n",
" y[:] = f(el, args...)\n",
" el.attributes[field] = copy(orig)\n",
" end\n",
" jac = ForwardDiff.forwarddiff_jacobian(helper!, Float64, fadtype=:dual, n=dim*nnodes, m=dim*nnodes)\n",
" return jac(el.attributes[field][:])\n",
" end\n",
" return jacobian\n",
"end\n",
"\n",
"function linearize!(f::Function, el::JuliaFEM.Element, field::ASCIIString, target::ASCIIString)\n",
" el.attributes[target][:] = 0.0\n",
" dim, nnodes = size(el.attributes[field])\n",
" function helper!(x, y)\n",
" orig = copy(el.attributes[field])\n",
" el.attributes[field] = reshape(x, dim, nnodes)\n",
" y[:] = f(el)\n",
" el.attributes[field] = copy(orig)\n",
" end\n",
" jac! = ForwardDiff.forwarddiff_jacobian!(helper!, Float64, fadtype=:dual, n=dim*nnodes, m=dim*nnodes)\n",
" jac!(el.attributes[field][:], el.attributes[target])\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"8x8 Array{Float64,2}:\n",
" 149.721 55.55 84.679 36.65 … -55.55 -136.278 -36.65 \n",
" 55.55 329.69 42.75 167.935 -172.941 -42.8 -324.684\n",
" 84.679 42.75 185.321 105.05 -123.05 -73.522 -24.75 \n",
" 36.65 167.935 105.05 340.54 -344.759 -24.8 -163.716\n",
" -98.122 -55.5 -196.478 -116.9 135.8 76.233 36.6 \n",
" -55.55 -172.941 -123.05 -344.759 … 352.922 42.8 164.778\n",
" -136.278 -42.8 -73.522 -24.8 42.8 133.567 24.8 \n",
" -36.65 -324.684 -24.75 -163.716 164.778 24.8 323.622"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"integrate(linearize(calc_residual_vector_integrand, \"displacement\"))(e)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"8x8 Array{Float64,2}:\n",
" 149.721 55.55 84.679 36.65 … -55.55 -136.278 -36.65 \n",
" 55.55 329.69 42.75 167.935 -172.941 -42.8 -324.684\n",
" 84.679 42.75 185.321 105.05 -123.05 -73.522 -24.75 \n",
" 36.65 167.935 105.05 340.54 -344.759 -24.8 -163.716\n",
" -98.122 -55.5 -196.478 -116.9 135.8 76.233 36.6 \n",
" -55.55 -172.941 -123.05 -344.759 … 352.922 42.8 164.778\n",
" -136.278 -42.8 -73.522 -24.8 42.8 133.567 24.8 \n",
" -36.65 -324.684 -24.75 -163.716 164.778 24.8 323.622"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"linearize(integrate(calc_residual_vector_integrand), \"displacement\")(e)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"8x8 Array{Float64,2}:\n",
" 149.721 55.55 84.679 36.65 … -55.55 -136.278 -36.65 \n",
" 55.55 329.69 42.75 167.935 -172.941 -42.8 -324.684\n",
" 84.679 42.75 185.321 105.05 -123.05 -73.522 -24.75 \n",
" 36.65 167.935 105.05 340.54 -344.759 -24.8 -163.716\n",
" -98.122 -55.5 -196.478 -116.9 135.8 76.233 36.6 \n",
" -55.55 -172.941 -123.05 -344.759 … 352.922 42.8 164.778\n",
" -136.278 -42.8 -73.522 -24.8 42.8 133.567 24.8 \n",
" -36.65 -324.684 -24.75 -163.716 164.778 24.8 323.622"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"linearize(integrate(calc_residual_vector_integrand), e, \"displacement\")"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"8x8 Array{Float64,2}:\n",
" 149.721 55.55 84.679 36.65 … -55.55 -136.278 -36.65 \n",
" 55.55 329.69 42.75 167.935 -172.941 -42.8 -324.684\n",
" 84.679 42.75 185.321 105.05 -123.05 -73.522 -24.75 \n",
" 36.65 167.935 105.05 340.54 -344.759 -24.8 -163.716\n",
" -98.122 -55.5 -196.478 -116.9 135.8 76.233 36.6 \n",
" -55.55 -172.941 -123.05 -344.759 … 352.922 42.8 164.778\n",
" -136.278 -42.8 -73.522 -24.8 42.8 133.567 24.8 \n",
" -36.65 -324.684 -24.75 -163.716 164.778 24.8 323.622"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"linearize!(integrate(calc_residual_vector_integrand), e, \"displacement\", \"displacement tangent stiffness\")\n",
"e.attributes[\"displacement tangent stiffness\"]"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"8x8 Array{Float64,2}:\n",
" 149.721 55.55 84.679 36.65 … -55.55 -136.278 -36.65 \n",
" 55.55 329.69 42.75 167.935 -172.941 -42.8 -324.684\n",
" 84.679 42.75 185.321 105.05 -123.05 -73.522 -24.75 \n",
" 36.65 167.935 105.05 340.54 -344.759 -24.8 -163.716\n",
" -98.122 -55.5 -196.478 -116.9 135.8 76.233 36.6 \n",
" -55.55 -172.941 -123.05 -344.759 … 352.922 42.8 164.778\n",
" -136.278 -42.8 -73.522 -24.8 42.8 133.567 24.8 \n",
" -36.65 -324.684 -24.75 -163.716 164.778 24.8 323.622"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"integrate!(linearize(calc_residual_vector_integrand, \"displacement\"), e, \"displacement tangent stiffness\")\n",
"e.attributes[\"displacement tangent stiffness\"]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Validations"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Converged in 6 iterations.\n"
]
},
{
"data": {
"text/plain": [
"2x4 Array{Float64,2}:\n",
" 0.0 -0.399145 -0.0722858 0.0\n",
" 0.0 -2.17799 -2.22224 0.0"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"free_dofs = [3, 4, 5, 6]\n",
"u = zeros(2, 4)\n",
"du = zeros(2, 4)\n",
"F = [0 0; 0 0; 0 -2; 0 0]'\n",
"for i=1:10\n",
" e.attributes[\"displacement\"] = u\n",
" K = linearize(integrate(calc_residual_vector_integrand), \"displacement\")(e)\n",
" R = integrate(calc_residual_vector_integrand)(e)\n",
" du[free_dofs] = K[free_dofs, free_dofs] \\ -(R - F)[free_dofs]\n",
" u += du\n",
" if norm(du) < 1.0e-9\n",
" println(\"Converged in $i iterations.\")\n",
" break\n",
" end\n",
"end\n",
"u # -2.222244754401764"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 0.4.0-dev",
"language": "julia",
"name": "julia-0.4"
},
"language_info": {
"name": "julia",
"version": "0.4.0"
}
},
"nbformat": 4,
"nbformat_minor": 0
}
notebooks/2015-08-25-interpolation-of-field-variables.ipynb
+362
View File
File diff suppressed because one or more lines are too long
src/JuliaFEM.jl
+10
View File
@@ -4,10 +4,20 @@ module JuliaFEM
VERSION < v"0.4-" && using Docile VERSION < v"0.4-" && using Docile
using Lexicon using Lexicon
using Logging
@Logging.configure(level=DEBUG)
Logging.info("loading types")
include("types.jl") # type definitions
Logging.info("loading elements")
include("elements.jl") # elements
include("math.jl") # basic mathematical operations
include("elasticity_solver.jl") include("elasticity_solver.jl")
include("xdmf.jl") include("xdmf.jl")
include("abaqus_reader.jl") include("abaqus_reader.jl")
include("interfaces.jl") include("interfaces.jl")
export set_coordinates, get_coordinates, set_material
end # module end # module
src/abaqus_reader.jl
+9 -17
View File
@@ -1,17 +1,12 @@
# This file is a part of JuliaFEM. # This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
module abaqus_reader eldims = Dict(
"C3D10" => 10,
"C3D4" => 4)
using Logging
@Logging.configure(level=DEBUG)
VERSION < v"0.4-" && using Docile
eldims = Dict({"C3D10" => 10})
global handlers = Dict() global handlers = Dict()
""" """
Register new handler for parser Register new handler for parser
""" """
@@ -29,13 +24,12 @@ end
function parse_header(header_line) function parse_header(header_line)
args = map(s -> strip(s), split(header_line, ",")) args = map(s -> strip(s), split(header_line, ","))
args[1] = strip(args[1], '*') args[1] = strip(args[1], '*')
d = Dict({"section" => args[1]}) d = Dict("section" => args[1], "options" => Dict())
options = Dict() options = d["options"]
for k in args[2:end] for k in args[2:end]
args2 = split(k, "=") args2 = split(k, "=")
options[args2[1]] = args2[2] options[args2[1]] = args2[2]
end end
d["options"] = options
return d return d
end end
@@ -56,7 +50,7 @@ function parse_element_section(model, header, data)
end end
eldim = eldims[eltype] eldim = eldims[eltype]
m = matchall(r"[0-9]+", data) m = matchall(r"[0-9]+", data)
m = map(integer, m) m = map((s) -> parse(Int, s), m)
elements = create_or_get(model, "elements") elements = create_or_get(model, "elements")
m = reshape(m, eldim+1, round(Int, length(m)/(eldim+1))) m = reshape(m, eldim+1, round(Int, length(m)/(eldim+1)))
nel = size(m)[2] nel = size(m)[2]
@@ -80,7 +74,7 @@ function parse_nodeset_section(model, header, data)
nset_name = header["options"]["NSET"] nset_name = header["options"]["NSET"]
Logging.debug("Creating node set $nset_name") Logging.debug("Creating node set $nset_name")
m = matchall(r"[0-9]+", data) m = matchall(r"[0-9]+", data)
node_ids = map(integer, m) node_ids = map((s) -> parse(Int, s), m)
nsets = create_or_get(model, "nsets") nsets = create_or_get(model, "nsets")
nsets[nset_name] = Int64[] nsets[nset_name] = Int64[]
for j in node_ids for j in node_ids
@@ -110,10 +104,10 @@ function parse_abaqus(fid)
end end
for line in eachline(fid) for line in eachline(fid)
if beginswith(line, "**") if startswith(line, "**")
continue continue
end end
if beginswith(line, "*") if startswith(line, "*")
process_section(section) process_section(section)
header = parse_header(line) header = parse_header(line)
Logging.debug("Found ", header["section"], " section") Logging.debug("Found ", header["section"], " section")
@@ -131,5 +125,3 @@ add_handler("NODE", parse_node_section)
add_handler("ELEMENT", parse_element_section) add_handler("ELEMENT", parse_element_section)
add_handler("NSET", parse_nodeset_section) add_handler("NSET", parse_nodeset_section)
end
src/elasticity_solver.jl
+85 -82
View File
@@ -3,6 +3,8 @@
module elasticity_solver module elasticity_solver
using ForwardDiff
using Logging using Logging
@Logging.configure(level=INFO) @Logging.configure(level=INFO)
@@ -12,103 +14,104 @@ VERSION < v"0.4-" && using Docile
# directly if needed or using general interface combining data model and # directly if needed or using general interface combining data model and
# solver. # solver.
""" """
Interpolate field variable using basis functions f for point ip. This is dummy function. Testing doctests and documentation.
This function tries to be as general as possible and allows interpolating
lot of different fields.
Parameters Parameters
---------- ----------
field :: Array{Number, dim} x : Array{Float64, 1}
Field variable
basis :: Function Returns
Basis functions -------
ip :: Array{Number, 1} Array{float64, 1}
Point to interpolate x + 1
Notes
-----
This is dummy function
Raises
------
Exception
if things are not going right
Examples
--------
>>> a = [1.0, 2.0, 3.0]
>>> dummy(a)
[2.0, 3.0, 4.0]
""" """
function interpolate{T<:Real}(field::Array{T,1}, basis::Function, ip) function dummy(a)
result = dot(field, basis(ip)) # not doing anything useful.
return result return a+1
end
function interpolate{T<:Real}(field::Array{T,2}, basis::Function, ip)
m, n = size(field)
bip = basis(ip)
tmp = size(bip)
if length(tmp) == 1
ndim = 1
nnodes = tmp[1]
else
ndim, nnodes = size(bip)
end
if ndim == 1
if n == nnodes
result = field * bip
elseif m == nnodes
result = field' * bip
end
else
if n == nnodes
result = bip' * field
elseif m == nnodes
result = bip' * field'
end
end
if length(result) == 1
result = result[1]
end
return result
end end
""" """
Calculate local tangent stiffness matrix and residual force vector R = T - F Calculate local tangent stiffness matrix and residual force vector
R = T - F for elasticity problem.
Parameters
----------
X : Element coordinates
u : Displacement field
R : Residual force vector
K : Tangent stiffness matrix
basis : Basis functions
dbasis : Derivative of basis functions
lambda : Material parameter
mu : Material parameter
ipoints : integration points
iweights : integration weights
Returns
-------
None
Notes
-----
If material parameters are given in list, they are interpolated to gauss
points using shape functions.
""" """
function calc_local_matrices!(X, u, R, Kt, N, dNdchi, lambda_, mu_, ipoints, iweights) function calc_local_matrices!(X, u, R, K, basis, dbasis, lambda_, mu_, ipoints, iweights)
dim, nnodes = size(X) dim, nnodes = size(X)
I = eye(dim) I = eye(dim)
R[:,:] = 0.0 R[:,:] = 0.0
Kt[:,:] = 0.0
dF = zeros(dim, dim) #dF = zeros(dim, dim)
for m = 1:length(iweights) function calc_R!(u, R)
w = iweights[m] for m = 1:length(iweights)
chi = ipoints[m, :] w = iweights[m]
# interpolate material parameters from element node fields xi = ipoints[m, :]
#lambda = (lambda_*N(chi))[1] # calculate material parameters
#mu = (mu_*N(chi))[1] lambda = typeof(lambda_) == Float64 ? lambda_ : dot(lambda_, basis(xi))
# Jt = X*dNdchi(chi) mu = typeof(mu_) == Float64 ? mu_ : dot(mu_, basis(xi))
#@debug("Jt:\n",Jt) Jt = X*dbasis(xi)
lambda = interpolate(lambda_, N, chi) detJ = det(Jt)
mu = interpolate(mu_, N, chi) dbasisdX = dbasis(xi)*inv(Jt)
Jt = interpolate(X, dNdchi, chi)
detJ = det(Jt)
deltaN = inv(Jt)*dNdchi(chi)'
delta_u = u*deltaN'
F = I + delta_u # Deformation gradient
E = 1/2*(delta_u' + delta_u + delta_u'*delta_u) # Green-Lagrange strain tensor
S = lambda*trace(E)*I + 2*mu*E # PK2 stress tensor
P = F*S # PK1 stress tensor
R[:,:] += w*P*deltaN*detJ
for p = 1:nnodes gradu = u*dbasisdX
for i = 1:dim F = I + gradu # Deformation gradient
dF[:,:] = 0.0 E = 1/2*(gradu' + gradu + gradu'*gradu) # Green-Lagrange strain tensor
dF[i,:] = deltaN[:,p] S = lambda*trace(E)*I + 2*mu*E # PK2 stress tensor
dE = 1/2*(F'*dF + dF'*F) P = F*S # PK1 stress tensor
dS = lambda*trace(dE)*I + 2*mu*dE
dP = dF*S + F*dS R[:,:] += w*P*dbasisdX'*detJ
for q = 1:nnodes
for j = 1:dim
Kt[dim*(p-1)+i,dim*(q-1)+j] += w*(dP[j,:]*deltaN[:,q])[1]*detJ
end
end
end
end end
end end
# herlper for tangent stiffness matrix
function R!(u, R)
R[:] = 0
calc_R!(reshape(u, dim, nnodes), reshape(R, dim, nnodes))
#calc_Wext!(reshape(u, 2, 4), reshape(R, 2, 4))
end
Jacobian = ForwardDiff.forwarddiff_jacobian(R!, Float64, fadtype=:dual, n=dim*nnodes, m=dim*nnodes)
K[:, :] = Jacobian(reshape(u, dim*nnodes))
R!(reshape(u, dim*nnodes), reshape(R, dim*nnodes))
end end
src/elements.jl
+235
View File
@@ -0,0 +1,235 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
abstract Element
"""
Get jacobian of element evaluated at point xi
"""
function get_jacobian(el::Element, xi)
dbasisdxi(xi) = get_dbasisdxi(el, xi)
X = get_coordinates(el)
J = interpolate(X, dbasisdxi, xi)'
return J
end
"""
Evaluate partial derivatives of basis function w.r.t
material description X, i.e. dbasis/dX
"""
function get_dbasisdX(el::Element, xi)
dbasisdxi = get_dbasisdxi(el, xi)
J = get_jacobian(el, xi)
dbasisdxi*inv(J)
end
"""
Return coordinates of element in array of size dim x nnodes
"""
function get_coordinates(el::Element)
el.coordinates
end
"""
Set coordinates for element
"""
function set_coordinates(el::Element, coordinates)
el.coordinates = coordinates
end
"""
Get element id
"""
function get_element_id(el::Element)
el.id
end
### Lagrange family ###
abstract CG <: Element # Lagrange element family
"""
Create new Lagrange element
FIXME: this is not working
LoadError: error compiling anonymous: type definition not allowed inside a local scope
It's the for loop which is causing problems. See
https://github.com/JuliaLang/julia/issues/10555
"""
function create_lagrange_element(element_name, X, P, dP)
@eval begin
nnodes, dim = size(X)
A = zeros(nnodes, nnodes)
for i=1:nnodes
A[i,:] = P(X[i,:])
end
invA = inv(A)'
type $element_name
element_id :: Int
node_ids :: Array{Int, 1}
coordinates :: Array{Float64, 2}
fields :: Dict{ASCIIString, Any}
end
function $element_name(element_id, node_ids)
coordinates = zeros(dim, nnodes)
fields = Dict{ASCIIString, Any}()
$element_name(element_id, node_ids, coordinates, fields)
end
function $element_name(element_id, node_ids, coordinates)
fields = Dict{ASCIIString, Any}()
$element_name(element_id, node_ids, coordinates, fields)
end
function get_basis(el::$element_name, xi)
invA*P(xi)
end
function get_dbasisdxi(el::$element_name, xi)
invA*dP(xi)
end
$element_name
end
end
# 0d Lagrange elements
"""
1 node point element
"""
type Point1 <: CG
element_id :: Int
node_ids :: Array{Int, 1}
coordinates :: Array{Float64, 2}
fields :: Dict{ASCIIString, Any}
end
# 1d Lagrange elements
"""
2 node linear line element
"""
type Seg2 <: CG
element_id :: Int
node_ids :: Array{Int, 1}
coordinates :: Array{Float64, 2}
fields :: Dict{ASCIIString, Any}
end
# X = [-1.0 1.0]'
# P = (xi) -> [1.0 xi[1]]'
# dP = (xi) -> [0.0 1.0]'
# create_lagrange_element(:Seg2, X, P, dP)
"""
3 node quadratic line element
"""
type Seg2 <: CG
element_id :: Int
node_ids :: Array{Int, 1}
coordinates :: Array{Float64, 2}
fields :: Dict{ASCIIString, Any}
end
#X = [-1.0 1.0 0.0]'
#P = (xi) -> [1.0 xi[1] xi[1]^2]'
#dP = (xi) -> [0.0 1.0 2*xi[1]]'
#create_lagrange_element(:Seg3, X, P, dP)
# 2d Lagrange elements
"""
4 node bilinear quadrangle element
"""
type Quad4 <: CG
element_id :: Int
node_ids :: Array{Int, 1}
coordinates :: Array{Float64, 2}
fields :: Dict{ASCIIString, Any}
end
function get_basis(el::Quad4, xi)
[(1-xi[1])*(1-xi[2])/4
(1+xi[1])*(1-xi[2])/4
(1+xi[1])*(1+xi[2])/4
(1-xi[1])*(1+xi[2])/4]
end
function get_dbasisdxi(el::Quad4, xi)
[-(1-xi[2])/4.0 -(1-xi[1])/4.0
(1-xi[2])/4.0 -(1+xi[1])/4.0
(1+xi[2])/4.0 (1+xi[1])/4.0
-(1+xi[2])/4.0 (1-xi[1])/4.0]
end
#X = [
# -1.0 -1.0
# 1.0 -1.0
# 1.0 1.0
# -1.0 1.0]
#P = (xi) -> [
# 1.0
# xi[1]
# xi[2]
# xi[1]*xi[2]]
#dP = (xi) -> [
# 0.0 0.0
# 1.0 0.0
# 0.0 1.0
# xi[2] xi[1]]
#create_lagrange_element(:Quad4, X, P, dP)
# 3d Lagrange elements
"""
10 node quadratic tethahedron
"""
type Tet10 <: CG
element_id :: Int
node_ids :: Array{Int, 1}
coordinates :: Array{Float64, 2}
fields :: Dict{ASCIIString, Any}
end
# X = [
# 0.0 0.0 0.0
# 1.0 0.0 0.0
# 0.0 1.0 0.0
# 0.0 0.0 1.0
# 0.5 0.0 0.0
# 0.5 0.5 0.0
# 0.0 0.5 0.0
# 0.0 0.0 0.5
# 0.5 0.0 0.5
# 0.0 0.5 0.5]
# P(xi) = [
# 1
# xi[1]
# xi[2]
# xi[3]
# xi[1]^2
# xi[2]^2
# xi[3]^2
# xi[1]*xi[2]
# xi[2]*xi[3]
# xi[3]*xi[1]]
# dP(xi) = [
# 0 0 0
# 1 0 0
# 0 1 0
# 0 0 1
# 2*xi[1] 0 0
# 0 2*xi[2] 0
# 0 0 2*xi[3]
# xi[2] xi[1] 0
# 0 xi[3] xi[2]
# xi[3] 0 xi[1]
# ]
#create_lagrange_element(:Tet10, X, P, dP)
src/math.jl
+186
View File
@@ -0,0 +1,186 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
"""
This module contains math stuff, including interpolation, integration, linearization, ...
"""
using ForwardDiff
export interpolate, integrate, linearize
"""
Interpolate field variable using basis functions f for point ip.
This function tries to be as general as possible and allows interpolating
lot of different fields.
Parameters
----------
field :: Array{Number, dim}
Field variable
basis :: Function
Basis functions
ip :: Array{Number, 1}
Point to interpolate
"""
function interpolate(field::Float64, basis::Function, ip::Array{Float64,1})
# dummy function, unable to interpolate scalar value!
return field
end
function interpolate{T<:Real}(field::Array{T,1}, basis::Function, ip)
result = dot(field, basis(ip))
return result
end
function interpolate{T<:Real}(field::Array{T,2}, basis::Function, ip)
m, n = size(field)
bip = basis(ip)
tmp = size(bip)
if length(tmp) == 1
ndim = 1
nnodes = tmp[1]
else
ndim, nnodes = size(bip)
end
if ndim == 1
if n == nnodes
result = field * bip
elseif m == nnodes
result = field' * bip
end
else
if n == nnodes
result = bip' * field
elseif m == nnodes
result = bip' * field'
end
end
if length(result) == 1
result = result[1]
end
return result
end
function interpolate(e::Element, field::ASCIIString, x::Array{Float64,1}; derivative=false)
basis = derivative ? get_dbasisdxi(e) : get_basis(e)
return interpolate(e.attributes[field], basis, x)
end
"""
Linearize function f w.r.t some given field, i.e. calculate dR/du
Parameters
----------
f::Function
(possibly) nonlinear function to linearize
field::ASCIIString
field variable
Returns
-------
Array{Float64, 2}
jacobian / "tangent stiffness matrix"
"""
function linearize(f::Function, el::Element, field::ASCIIString)
dim, nnodes = size(el.attributes[field])
function helper!(x, y)
orig = copy(el.attributes[field])
el.attributes[field] = reshape(x, dim, nnodes)
y[:] = f(el)
el.attributes[field] = copy(orig)
end
jac = ForwardDiff.forwarddiff_jacobian(helper!, Float64, fadtype=:dual, n=dim*nnodes, m=dim*nnodes)
return jac(el.attributes[field][:])
end
"""
This version returns another function which can be then evaluated against field
"""
function linearize(f::Function, field::ASCIIString)
function jacobian(el::Element, args...)
dim, nnodes = size(el.attributes[field])
function helper!(x, y)
orig = copy(el.attributes[field])
el.attributes[field] = reshape(x, dim, nnodes)
y[:] = f(el, args...)
el.attributes[field] = copy(orig)
end
jac = ForwardDiff.forwarddiff_jacobian(helper!, Float64, fadtype=:dual, n=dim*nnodes, m=dim*nnodes)
return jac(el.attributes[field][:])
end
return jacobian
end
"""
In-place version, no additional garbage collection.
"""
function linearize!(f::Function, el::Element, field::ASCIIString, target::ASCIIString)
el.attributes[target][:] = 0.0
dim, nnodes = size(el.attributes[field])
function helper!(x, y)
orig = copy(el.attributes[field])
el.attributes[field] = reshape(x, dim, nnodes)
y[:] = f(el)
el.attributes[field] = copy(orig)
end
jac! = ForwardDiff.forwarddiff_jacobian!(helper!, Float64, fadtype=:dual, n=dim*nnodes, m=dim*nnodes)
jac!(el.attributes[field][:], el.attributes[target])
end
"""
Integrate f over element using Gaussian quadrature rules.
Parameters
----------
el::Element
well defined element
f::Function
Function to integrate
"""
function integrate(f::Function, el::Element)
target = []
for ip in el.integration_points
J = interpolate(el, "coordinates", ip.xi; derivative=true)
push!(target, ip.weight*f(el, ip)*det(J))
end
return sum(target)
end
#function integrate(f::Function, integration_points::Array{IntegrationPoint, 1}, Xargs...)
# target = []
# for ip in integration_points
# J = interpolate(el, "coordinates", ip.xi; derivative=true)
# push!(target, ip.weight*f(ip, args...)*det(J))
# end
# return sum(target)
#end
"""
This version returns a function which must be operated with element e
"""
function integrate(f::Function)
function integrate(el::Element)
target = []
for ip in el.integration_points
J = interpolate(el, "coordinates", ip.xi; derivative=true)
push!(target, ip.weight*f(el, ip)*det(J))
end
return sum(target)
end
return integrate
end
"""
This version saves results inplace to target, garbage collection free
"""
function integrate!(f::Function, el::Element, target)
# set target to zero
el.attributes[target][:] = 0.0
for ip in el.integration_points
J = interpolate(el, "coordinates", ip.xi; derivative=true)
el.attributes[target][:,:] += ip.weight*f(el, ip)*det(J)
end
end
src/types.jl
+52
View File
@@ -0,0 +1,52 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
export IntegrationPoint, Element, Assembly, FunctionSpace
"""
Integration point
xi :: Array{Float64, 1}
(dimensionless) coordinates of integration point
weight :: Float64
Integration weight
attributes :: Dict{ASCIIString, Any}
This is used to save internal variables of IP needed e.g. for incremental
material models.
"""
type IntegrationPoint
xi :: Array{Float64, 1}
weight :: Float64
attributes :: Dict{ASCIIString, Any}
end
type FunctionSpace
basis :: Function
dbasis :: Function
end
abstract Element
#type Element
# id :: Int
# node_ids :: Array{Int, 1}
# shape_functions :: FunctionSpace
# integration_points :: Array{IntegrationPoint, 1}
# attributes :: Dict{ASCIIString, Any}
#end
type Assembly
# LHS
I :: Array{Int64, 1}
J :: Array{Int64, 1}
A :: Array{Float64, 1}
# RHS
i :: Array{Int64, 1}
b :: Array{Float64, 1}
# global dofs for each element
gdofs :: Dict{Int64, Array{Int64, 1}}
end
Assembly() = Assembly(Int64[], Int64[], Float64[], Int64[], Float64[], Dict{Int64,Array{Int64,1}}())
Assembly(gdofs::Dict{Int64,Array{Int64,1}}) = Assembly(Int64[], Int64[], Float64[], Int64[], Float64[], gdofs)
src/xdmf.jl
+51 -22
View File
@@ -1,17 +1,8 @@
# This file is a part of JuliaFEM. # This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
module xdmf
using Logging
@Logging.configure(level=INFO)
using LightXML using LightXML
VERSION < v"0.4-" && using Docile
# i add docstrings later
# element codes: http://www.paraview.org/pipermail/paraview/2013-July/028859.html # element codes: http://www.paraview.org/pipermail/paraview/2013-July/028859.html
# > from ./VTK/ThirdParty/xdmf2/vtkxdmf2/libsrc/XdmfTopology.h # > from ./VTK/ThirdParty/xdmf2/vtkxdmf2/libsrc/XdmfTopology.h
# > # >
@@ -84,29 +75,68 @@ function xdmf_new_grid(temporal_collection; time=0)
return grid return grid
end end
#function xdmf_new_mesh(grid, X, elmap)
# geometry = new_child(grid, "Geometry")
# set_attribute(geometry, "Type", "XYZ")
# dataitem = new_child(geometry, "DataItem")
# set_attribute(dataitem, "DataType", "Float")
# set_attribute(dataitem, "Dimensions", length(X))
# set_attribute(dataitem, "Format", "XML")
# set_attribute(dataitem, "Precision", "4")
# add_text(dataitem, join(X, " "))
# topology = new_child(grid, "Topology")
# set_attribute(topology, "Dimensions", "1")
# set_attribute(topology, "Type", "Mixed")
# dataitem = new_child(topology, "DataItem")
# set_attribute(dataitem, "DataType", "Int")
# set_attribute(dataitem, "Dimensions", length(elmap))
# set_attribute(dataitem, "Format", "XML")
# set_attribute(dataitem, "Precision", 4)
# elmap2 = copy(elmap)
# elmap2[2:end,:] -= 1
# add_text(dataitem, join(elmap2, " "))
#end
function xdmf_new_mesh(grid, X, elmap) function xdmf_new_mesh(grid, X, elmap)
dim, nnodes = size(X)
geometry = new_child(grid, "Geometry") geometry = new_child(grid, "Geometry")
set_attribute(geometry, "Type", "XYZ") set_attribute(geometry, "Type", "XYZ")
dataitem = new_child(geometry, "DataItem") dataitem = new_child(geometry, "DataItem")
set_attribute(dataitem, "DataType", "Float") set_attribute(dataitem, "DataType", "Float")
set_attribute(dataitem, "Dimensions", length(X)) set_attribute(dataitem, "Dimensions", "$nnodes $dim")
set_attribute(dataitem, "Format", "XML") set_attribute(dataitem, "Format", "XML")
set_attribute(dataitem, "Precision", "4") set_attribute(dataitem, "Precision", 8)
add_text(dataitem, join(X, " ")) #add_text(dataitem, join(X, " "))
s = "\n"
for i=1:nnodes
s *= "\t\t" * join(X[:,i], " ") * "\n"
end
s *= " "
add_text(dataitem, s)
topology = new_child(grid, "Topology")
set_attribute(topology, "Dimensions", "1")
set_attribute(topology, "Type", "Mixed")
dataitem = new_child(topology, "DataItem")
set_attribute(dataitem, "DataType", "Int")
set_attribute(dataitem, "Dimensions", length(elmap))
set_attribute(dataitem, "Format", "XML")
set_attribute(dataitem, "Precision", 4)
elmap2 = copy(elmap) elmap2 = copy(elmap)
elmap2[2:end,:] -= 1 elmap2[2:end,:] -= 1
add_text(dataitem, join(elmap2, " ")) dim, nelements = size(elmap2)
topology = new_child(grid, "Topology")
#set_attribute(topology, "Dimensions", "1")
set_attribute(topology, "TopologyType", "Mixed")
set_attribute(topology, "NumberOfElements", nelements)
dataitem = new_child(topology, "DataItem")
set_attribute(dataitem, "DataType", "Int")
set_attribute(dataitem, "Dimensions", "$nelements $dim")
set_attribute(dataitem, "Format", "XML")
set_attribute(dataitem, "Precision", 8)
s = "\n"
for i=1:nelements
s *= "\t\t" * join(elmap2[:,i], " ") * "\n"
end
add_text(dataitem, s)
#add_text(dataitem, join(elmap2, " "))
end end
function xdmf_new_field(grid, name, source, data) function xdmf_new_field(grid, name, source, data)
loc = Dict("elements" => "Cell", loc = Dict("elements" => "Cell",
"nodes" => "Node") "nodes" => "Node")
@@ -149,4 +179,3 @@ function xdmf_save_model(xdoc, filename)
save_file(xdoc, filename) save_file(xdoc, filename)
end end
end
test/test_abaqus_reader.jl
+45 -29
View File
@@ -5,45 +5,61 @@ using FactCheck
using Logging using Logging
@Logging.configure(level=INFO) @Logging.configure(level=INFO)
using JuliaFEM.abaqus_reader: parse_abaqus, parse_element_section #using JuliaFEM.abaqus_reader: parse_abaqus, parse_element_section
include(Pkg.dir("JuliaFEM")*"/src/abaqus_reader.jl")
facts("test import abaqus model") do facts("test import abaqus model") do
# FIXME: get_test_data() # FIXME: get_test_data()
fid = open(Pkg.dir("JuliaFEM")*"/geometry/3d_beam/palkki.inp") fid = open(Pkg.dir("JuliaFEM")*"/geometry/3d_beam/palkki.inp")
model = parse_abaqus(fid) model = parse_abaqus(fid)
close(fid) close(fid)
@fact length(model["nodes"]) => 298 @fact length(model["nodes"]) --> 298
@fact length(model["elements"]) => 120 @fact length(model["elements"]) --> 120
@fact length(model["elsets"]["Body1"]) => 120 @fact length(model["elsets"]["Body1"]) --> 120
@fact length(model["nsets"]["SUPPORT"]) => 9 @fact length(model["nsets"]["SUPPORT"]) --> 9
@fact length(model["nsets"]["LOAD"]) => 9 @fact length(model["nsets"]["LOAD"]) --> 9
@fact length(model["nsets"]["TOP"]) => 83 @fact length(model["nsets"]["TOP"]) --> 83
end end
facts("test that reader throws error when dimension information of elemenet is missing") do facts("test that reader throws error when dimension information of elemenet is missing") do
# *ELEMENT, TYPE=neverseenbefore, ELSET=Body1 # *ELEMENT, TYPE=neverseenbefore, ELSET=Body1
data = """ data = """
1, 243, 240, 191, 117, 245, 242, 244, 1, 243, 240, 191, 117, 245, 242, 244,
1, 2, 196 1, 2, 196
""" """
model = Dict() model = Dict()
header = Dict("section"=>"ELEMENT", "options" => Dict("TYPE" => "neverseenbefore", "ELSET"=>"Body1")) header = Dict("section"=>"ELEMENT", "options" => Dict("TYPE" => "neverseenbefore", "ELSET"=>"Body1"))
@fact_throws parse_element_section(model, header, data) @fact_throws parse_element_section(model, header, data)
end
facts("read element section") do
data = """
1, 243, 240, 191, 117, 245, 242, 244,
1, 2, 196
2, 204, 199, 175, 130, 207, 208, 209,
3, 4, 176
"""
model = Dict()
header = Dict("section" => "ELEMENT", "options" => Dict("TYPE" => "C3D10", "ELSET" => "BEAM"))
parse_element_section(model, header, data)
@fact length(model["elements"]) --> 2
@fact model["elements"][1] --> [243, 240, 191, 117, 245, 242, 244, 1, 2, 196]
@fact model["elements"][2] --> [204, 199, 175, 130, 207, 208, 209, 3, 4, 176]
end end
facts("test unknown handler warning message") do facts("test unknown handler warning message") do
fn = tempname() fn = tempname()
fid = open(fn, "w") fid = open(fn, "w")
testdata = """ testdata = """
*ELEMENT2, TYPE=C3D10, ELSET=Body1 *ELEMENT2, TYPE=C3D10, ELSET=Body1
1, 243, 240, 191, 117, 245, 242, 244, 1, 243, 240, 191, 117, 245, 242, 244,
1, 2, 196 1, 2, 196
""" """
write(fid, testdata) write(fid, testdata)
close(fid) close(fid)
fid = open(fn) fid = open(fn)
model = parse_abaqus(fid) model = parse_abaqus(fid)
close(fid) close(fid)
# empty model expected, parser doesn't know what to do with unknown section # empty model expected, parser doesn't know what to do with unknown section
@fact length(model) => 0 @fact length(model) --> 0
end end
test/test_elasticity_solver.jl
+119 -29
View File
@@ -140,35 +140,6 @@ facts("test solve elasticity increment, two elements") do
end end
using JuliaFEM.elasticity_solver: interpolate
facts("test interpolation of different field variables") do
N(xi) = [
(1-xi[1])*(1-xi[2])/4
(1+xi[1])*(1-xi[2])/4
(1+xi[1])*(1+xi[2])/4
(1-xi[1])*(1+xi[2])/4
]
dNdξ(ξ) = [-(1-ξ[2])/4.0 -(1-ξ[1])/4.0
(1-ξ[2])/4.0 -(1+ξ[1])/4.0
(1+ξ[2])/4.0 (1+ξ[1])/4.0
-(1+ξ[2])/4.0 (1-ξ[1])/4.0]
F1 = [36.0, 36.0, 36.0, 36.0]
F2 = [36.0 36.0 36.0 36.0]
F3 = F2'
F4 = [0.0 0.0; 10.0 0.0; 10.0 1.0; 0.0 1.0]'
F5 = F4'
F6 = [36, 36, 36, 36]
@fact interpolate(F1, N, [0.0, 0.0]) => 36.0
@fact interpolate(F2, N, [0.0, 0.0]) => 36.0
@fact interpolate(F3, N, [0.0, 0.0]) => 36.0
@fact interpolate(F4, N, [0.0, 0.0]) => [5.0; 0.5]
@fact interpolate(F5, N, [0.0, 0.0]) => [5.0; 0.5]
@fact interpolate(F5, dNdξ, [0.0, 0.0]) => [5.0 0.0; 0.0 0.5]
@fact interpolate(F6, N, [0.0, 0.0]) => 36
end
using JuliaFEM.elasticity_solver: assemble! using JuliaFEM.elasticity_solver: assemble!
@@ -280,3 +251,122 @@ facts("test that elimination of non-homogeneous dirichlet boundary conditions ra
I, J, V = findnz(A) I, J, V = findnz(A)
@fact_throws I, V = eliminate_boundary_conditions(dirichletbc, I, V) @fact_throws I, V = eliminate_boundary_conditions(dirichletbc, I, V)
end end
module TestElasticitySolver
using JuliaFEM.elasticity_solver: calc_local_matrices
facts("test solve one element model") do
X = [0.0 0.0; 10.0 0.0; 10.0 1.0; 0.0 1.0]'
F = [0 0; 0 0; 0 -2; 0 0]'
# Material properties
E = 90
nu = 0.25
mu = E/(2*(1+nu))
la = E*nu/((1+nu)*(1-2*nu))
la = 2*la*mu/(la + 2*mu)
u = zeros(2, 4)
du = zeros(2, 4)
R = zeros(2, 4)
K = zeros(8, 8)
basis(xi) = [
(1-xi[1])*(1-xi[2])/4
(1+xi[1])*(1-xi[2])/4
(1+xi[1])*(1+xi[2])/4
(1-xi[1])*(1+xi[2])/4]
dbasis(xi) = [-(1-xi[2])/4.0 -(1-xi[1])/4.0
(1-xi[2])/4.0 -(1+xi[1])/4.0
(1+xi[2])/4.0 (1+xi[1])/4.0
-(1+xi[2])/4.0 (1-xi[1])/4.0]
ipoints = 1/sqrt(3)*[-1 -1; 1 -1; 1 1; -1 1]
iweights = [1, 1, 1, 1]
free_dofs = [3, 4, 5, 6]
for i=1:10
calc_local_matrices!(X, u, R, K, basis, dbasis, la, mu, ipoints, iweights)
du[free_dofs] = K[free_dofs, free_dofs] \ -(R - F)[free_dofs]
u += du
if norm(du) < 1.0e-9
Logging.debug("Converged in $i iterations.")
break
end
end
# Tested against Elmer solution
Logging.debug("solution vector: \n $u")
@fact u[2, 3] --> roughly(-2.222244754401764)
norm1 = norm(u)
Logging.debug("norm of u: $(norm(u))")
# We rotate model a bit and make sure that L2 norm is same
phi = 30/180*pi
rmat = [
cos(phi) -sin(phi)
sin(phi) cos(phi)]
X = rmat*X
F = rmat*F
u = zeros(2, 4)
for i=1:10
calc_local_matrices!(X, u, R, K, basis, dbasis, la, mu, ipoints, iweights)
du[free_dofs] = K[free_dofs, free_dofs] \ -(R - F)[free_dofs]
u += du
if norm(du) < 1.0e-9
Logging.debug("Converged in $i iterations.")
break
end
end
Logging.debug("solution vector: \n $u")
Logging.debug("norm of u: $(norm(u))")
@fact norm(u) --> roughly(norm1)
# test two element model
X = [0.0 0.0; 5.0 0.0; 5.0 1.0; 0.0 1.0]'
u = zeros(2, 6)
du = zeros(2, 6)
R = zeros(2, 4)
K = zeros(8, 8)
ass1 = [9, 10, 1, 2, 5, 6, 11, 12]
ass2 = [1, 2, 3, 4, 7, 8, 5, 6]
free_dofs = collect(1:8)
F = [0 0; 0 0; 0 0; 0 -0.1; 0 0; 0 0]'
A = zeros(12, 12)
b = zeros(2, 6)
for i=1:1
Logging.debug("Iteration $i")
A[:,:] = 0.0
b[:] = 0.0
#Logging.debug("Assembling")
for ass in (ass1, ass2)
#Logging.debug("ass = $ass, u[ass] = $(u[ass])")
calc_local_matrices!(X, u[ass], R, K, basis, dbasis, la, mu, ipoints, iweights)
A[ass,ass] += K
b[ass] += R[:]
end
dump(round(A, 2))
println("K norm = $(norm(A[free_dofs, free_dofs]))")
du[free_dofs] = A[free_dofs, free_dofs] \ -(b - F)[free_dofs]
println("du = $du")
u += du
Logging.debug("Norm of du: $(norm(du))")
for ass in (ass1, ass2)
Logging.debug("Element displacement: $(reshape(u[ass], 2, 4))")
end
if norm(du) < 1.0e-9
Logging.debug("Converged in $i iterations.")
break
end
end
Logging.debug("solution vector: \n $u")
Logging.debug("norm of u: $(norm(u))")
@pending norm(u) --> :something
end
exitstatus()
end
test/test_math.jl
+33
View File
@@ -0,0 +1,33 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
using JuliaFEM: interpolate
using FactCheck
facts("test interpolation of different field variables") do
N(xi) = [
(1-xi[1])*(1-xi[2])/4
(1+xi[1])*(1-xi[2])/4
(1+xi[1])*(1+xi[2])/4
(1-xi[1])*(1+xi[2])/4
]
dNdξ(ξ) = [-(1-ξ[2])/4.0 -(1-ξ[1])/4.0
(1-ξ[2])/4.0 -(1+ξ[1])/4.0
(1+ξ[2])/4.0 (1+ξ[1])/4.0
-(1+ξ[2])/4.0 (1-ξ[1])/4.0]
F1 = [36.0, 36.0, 36.0, 36.0]
F2 = [36.0 36.0 36.0 36.0]
F3 = F2'
F4 = [0.0 0.0; 10.0 0.0; 10.0 1.0; 0.0 1.0]'
F5 = F4'
F6 = [36, 36, 36, 36]
@fact interpolate(F1, N, [0.0, 0.0]) --> 36.0
@fact interpolate(F2, N, [0.0, 0.0]) --> 36.0
@fact interpolate(F3, N, [0.0, 0.0]) --> 36.0
@fact interpolate(F4, N, [0.0, 0.0]) --> [5.0; 0.5]
@fact interpolate(F5, N, [0.0, 0.0]) --> [5.0; 0.5]
@fact interpolate(F5, dNdξ, [0.0, 0.0]) --> [5.0 0.0; 0.0 0.5]
@fact interpolate(F6, N, [0.0, 0.0]) --> 36
end