Newer
Older
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
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
program BallGeometry
implicit none
integer level,i,j,nrad,nz,nradp,ip
integer i1,i2,i3,i4,ij,k,np
double precision x1,x2,x3,y1,y2,y3,z1,z2,z3
double precision xne,yne,zne,xyzn
integer nnode,nnode0,nxt,nxb,nnodemax,nelem,nelemmax,nnodep
double precision xc,yc,zc,rad,radp,phi,theta
double precision x,y,z,xsp,dist,pi
double precision, dimension(:),allocatable::xp,yp,zp,xn,yn,zn
integer,dimension(:,:),allocatable::icon
pi=atan(1.)*4.
! xc yc zc is the centre of the ball
! rad is radius
xc=.501d0
yc=.501d0
zc=.495d0
xc=.5d0
yc=.5d0
! zc=.9d0
rad=.025d0
! reference size
rad=.05d0
level=9
nz=int(2.d0**level*rad*pi)+1
if ((nz/2)*2.eq.nz) nz=nz+1
nnodemax=nz*nz*2
nelemmax=nnodemax*3
allocate (xp(nnodemax),yp(nnodemax),zp(nnodemax))
allocate (xn(nnodemax),yn(nnodemax),zn(nnodemax))
allocate (icon(3,nelemmax))
nnode=0
nelem=0
ip=1
do j=1,nz
theta=-pi/2.+pi*float(j-1)/float(nz-1)
radp=rad*cos(theta)
z=rad*sin(theta)
nrad=1+int(radp*pi*2.d0*2.d0**level)
print*,nz,nrad
do i=1,nrad
phi=2.d0*pi*float(i-1)/float(nrad)
x=radp*sin(phi)
y=radp*cos(phi)
call add_point (nnode,x,y,z,xp,yp,zp,nnodemax)
! bottom half
if (j.ne.1.and.j.le.nz/2+1) then
i1=nnode
i2=nnode+1
if (i.eq.nrad) i2=nnodep+nradp+1
i3=nnodep+int((i*nradp)/nrad)+1
ip=nnodep+int(((i-1)*nradp)/nrad)+1
if (i3.ge.nnodep+nradp+1) i3=nnodep+1
call add_element (nelem,icon,i1,i2,i3,nelemmax)
if (i3.ne.ip) call add_element (nelem,icon,i1,i3,ip,nelemmax)
ip=i3
! top half (minus very top)
elseif (j.ne.1.and.j.ne.nz) then
i1=nnode
i2=nnode+1
if (i.eq.nrad) i2=nnode-nrad+1
i3=nnodep+int((i*nradp)/nrad)+1
if (i.eq.1) i3=nnodep+2
np=1
if (i.ne.1) np=i3-ip
do k=1,np
call add_element (nelem,icon,i3-k,i1,i3-k+1,nelemmax)
enddo
if (i.eq.nrad) call add_element (nelem,icon,i2,nnodep+1,nnode-nrad,nelemmax)
call add_element (nelem,icon,i1,i2,i3,nelemmax)
ip=i3
! very top
elseif (j.eq.nz) then
do k=1,nradp
i1=nnodep+k
i2=i1+1
if (i2.eq.nnode) i2=nnodep+1
i3=nnode
call add_element (nelem,icon,i1,i3,i2,nelemmax)
enddo
endif
enddo
nradp=nrad
nnodep=nnode-nrad
enddo
print*,nnode,nelem
!add normals
112 continue
do i=1,nelem
i1=icon(1,i)
i2=icon(2,i)
i3=icon(3,i)
if (i1.eq.i2 .or. i2.eq.i3 .or. i3.eq.i1)then
print*,'element ',i,' collapsed: ',icon(:,i)
do j=i+1,nelem
icon(:,j-1)=icon(:,j)
enddo
nelem=nelem-1
goto 112
endif
x1=xp(i1)
x2=xp(i2)
x3=xp(i3)
y1=yp(i1)
y2=yp(i2)
y3=yp(i3)
z1=zp(i1)
z2=zp(i2)
z3=zp(i3)
xne=(y2-y1)*(z3-z1)-(y3-y1)*(z2-z1)
yne=(z2-z1)*(x3-x1)-(z3-z1)*(x2-x1)
zne=(x2-x1)*(y3-y1)-(x3-x1)*(y2-y1)
xyzn=sqrt(xne**2+yne**2+zne**2)
xne=xne/xyzn
yne=yne/xyzn
zne=zne/xyzn
do j=1,3
ij=icon(j,i)
xn(ij)=xn(ij)+xne
yn(ij)=yn(ij)+yne
zn(ij)=zn(ij)+zne
enddo
enddo
do i=1,nnode
xyzn=sqrt(xn(i)**2+yn(i)**2+zn(i)**2)
xn(i)=xn(i)/xyzn
yn(i)=yn(i)/xyzn
zn(i)=zn(i)/xyzn
enddo
xp=xp+xc
yp=yp+yc
zp=zp+zc
open (71,file='surface-1.txt',status='unknown')
write (71,*) nnode,nelem
do i=1,nnode
write (71,*) xp(i),yp(i),zp(i),xn(i),yn(i),zn(i)
enddo
do i=1,nelem
write (71,*) (icon(j,i),j=1,3)
enddo
close (71)
open(unit=30,file='Mesh.vtk')
write(30,'(a)')'# vtk DataFile Version 3.0'
write(30,'(a)')'surface'
write(30,'(a)')'ASCII'
write(30,'(a)')'DATASET UNSTRUCTURED_GRID'
write(30,'(a7,i10,a6)')'POINTS ',nnode,' float'
do i=1,nnode
write(30,'(3f16.11)') xp(i),yp(i),zp(i)
enddo
write(30,'(A6, 2I10)') 'CELLS ',nelem,4*nelem
do i=1,nelem
write(30,'(9I10)')3,icon(1:3,i)-1
enddo
write(30,'(A11, I10)') 'CELL_TYPES ',nelem
do i=1,nelem
write(30,'(I2)')5 ! octree (8 nodes)
enddo
write(30,'(a11,i10)')'POINT_DATA ',nnode
write(30,'(a)')'VECTORS norm float'
do i=1,nnode
write(30,'(3f18.13)')xn(i),yn(i),zn(i)
enddo
close(30)
deallocate (xp,yp,zp,icon)
end
!--------
subroutine add_point (nnode,x,y,z,xp,yp,zp,nnodemax)
double precision xp(nnodemax),yp(nnodemax),zp(nnodemax)
double precision x,y,z
nnode=nnode+1
if (nnode.gt.nnodemax) then
print*,nnode,nnodemax
stop 'needs to increase nnodemax'
endif
xp(nnode)=x
yp(nnode)=y
zp(nnode)=z
return
end
!----------
subroutine add_element (nelem,icon,i1,i2,i3,nelemmax)
integer nelem,nelemmax,icon(3,nelemmax),i1,i2,i3
nelem=nelem+1
if (nelem.gt.nelemmax) then
print*,nelem,nelemmax
stop 'nelemmax needs to be increased'
endif
icon(1,nelem)=i1
icon(2,nelem)=i2
icon(3,nelem)=i3
return
end