1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
|
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Python GetFEM++ interface
#
# Copyright (C) 2009 Yves Renard, Luis Saavedra.
#
# This file is a part of GetFEM++
#
# GetFEM++ is free software; you can redistribute it and/or modify it
# under the terms of the GNU Lesser General Public License as published
# by the Free Software Foundation; either version 2.1 of the License, or
# (at your option) any later version.
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
# or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
# License for more details.
# You should have received a copy of the GNU Lesser General Public License
# along with this program; if not, write to the Free Software Foundation,
# Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
#
############################################################################
""" test levelset.
This program is used to check that python-getfem is working. This is
also a good example of use of python-getfem..
$Id: check_levelset.py 4925 2015-03-30 16:12:19Z renard $
"""
import getfem as gf
import numpy as np
from scipy import rand
eps = 1.0/10
m = gf.Mesh('regular_simplices', np.arange(-1,1+eps,eps), np.arange(-1,1+eps,eps), 'degree', 2, 'noised')
#m = gf.Mesh('cartesian', np.arange(-1,1+eps,eps), np.arange(-1,1+eps,eps))
ls1 = gf.LevelSet(m, 2, 'y', 'x')
#ls1 = gf.LevelSet(m, 2, 'sqr(x) + sqr(y) - sqr(0.7)', 'x-.4')
#ls1 = gf.LevelSet(m, 2, 'x + y - 0.2') #, 'x-5')
#ls1 = gf.LevelSet(m, 2, 'x + y - 0.2', 'x-5')
#ls2 = gf.LevelSet(m, 2, '0.6*sqr(x) + sqr(y-0.1) - sqr(0.6)');
#ls3 = gf.LevelSet(m, 4, 'sqr(x) + sqr(y+.08) - sqr(0.05)');
ls2 = gf.LevelSet(m, 2, 'y+0.1', 'x')
ls3 = gf.LevelSet(m, 2, 'y-0.1', 'x')
mls = gf.MeshLevelSet(m)
mls.add(ls1)
if True:
mls.sup(ls1)
mls.add(ls1)
mls.add(ls2)
mls.add(ls2)
mls.add(ls2)
mls.add(ls3)
mls.adapt()
#print mls.linked_mesh()
lls = mls.levelsets()
cm = mls.cut_mesh()
ctip = mls.crack_tip_convexes()
mf = gf.MeshFem(m)
mf.set_classical_fem(1)
mfls = gf.MeshFem('levelset',mls,mf)
gf.memstats()
nbd = mfls.nbdof()
if True:
sl = gf.Slice(('none',), mls, 2);
U = rand(1,nbd);
sl.export_to_pos('slU.pos',mfls,U,'U')
mfls.export_to_pos('U.pos',U,'U')
cm.export_to_pos('cm.pos')
m.export_to_pos('m.pos')
else:
sl = gf.Slice(('none',), mls, 1);
for i in xrange(nbd):
U = np.zeros(nbd)
U[i] = 1
sl.export_to_pos('slU'+str(i)+'.pos',mfls,U,'U'+str(i))
mfls.export_to_pos('U'+str(i)+'.pos',U,'U'+str(i))
cm.export_to_pos('cm.pos')
m.export_to_pos('m.pos')
|
