Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
get_external_loads.f90
Go to the documentation of this file.
1 !*********************************************************************
2 !* Illinois Open Source License *
3 !* *
4 !* University of Illinois/NCSA *
5 !* Open Source License *
6 !* *
7 !* Copyright@2008, University of Illinois. All rights reserved. *
8 !* *
9 !* Developed by: *
10 !* *
11 !* Center for Simulation of Advanced Rockets *
12 !* *
13 !* University of Illinois *
14 !* *
15 !* www.csar.uiuc.edu *
16 !* *
17 !* Permission is hereby granted, free of charge, to any person *
18 !* obtaining a copy of this software and associated documentation *
19 !* files (the "Software"), to deal with the Software without *
20 !* restriction, including without limitation the rights to use, *
21 !* copy, modify, merge, publish, distribute, sublicense, and/or *
22 !* sell copies of the Software, and to permit persons to whom the *
23 !* Software is furnished to do so, subject to the following *
24 !* conditions: *
25 !* *
26 !* *
27 !* @ Redistributions of source code must retain the above copyright *
28 !* notice, this list of conditions and the following disclaimers. *
29 !* *
30 !* @ Redistributions in binary form must reproduce the above *
31 !* copyright notice, this list of conditions and the following *
32 !* disclaimers in the documentation and/or other materials *
33 !* provided with the distribution. *
34 !* *
35 !* @ Neither the names of the Center for Simulation of Advanced *
36 !* Rockets, the University of Illinois, nor the names of its *
37 !* contributors may be used to endorse or promote products derived *
38 !* from this Software without specific prior written permission. *
39 !* *
40 !* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, *
41 !* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES *
42 !* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND *
43 !* NONINFRINGEMENT. IN NO EVENT SHALL THE CONTRIBUTORS OR *
44 !* COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
45 !* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, *
46 !* ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
47 !* USE OR OTHER DEALINGS WITH THE SOFTWARE. *
48 !*********************************************************************
49 !* Please acknowledge The University of Illinois Center for *
50 !* Simulation of Advanced Rockets in works and publications *
51 !* resulting from this software or its derivatives. *
52 !*********************************************************************
53 SUBROUTINE get_external_loads(coords,jdltyp, &
54  adlmag,mdload,ndload,rex,nodes,numnp,&
55  mapselvolel,numelvol)
56 !
57 ! *----------------------------------------------GET_EXTERNAL_LOADS----*
58 ! | |
59 ! | Add surface tractions and Add body forces |
60 ! | |
61 ! | Notes (for details, see report): |
62 ! | |
63 ! | 1- Non-Uniform surface fluxes are supported |
64 ! | |
65 ! | 2- The face that the load is applied is designated by |
66 ! | load type in the |DLOAD card. For example : |
67 ! | *DLOAD, OP=NEW |
68 ! | SET_NAME, U1, 1.0 |
69 ! | describes a uniform pressure of magnitude 1.0 |
70 ! | to be applied to the face 1 (number taken from U1) |
71 ! | of all the elements in the set SET_NAME. |
72 ! | |
73 ! | 3- Load types mean, |
74 ! | U1, pressure onto the 1-2-3-4 face |
75 ! | U2, pressure onto the 5-8-7-6 face |
76 ! | U3, pressure onto the 1-5-6-2 face |
77 ! | U4, pressure onto the 2-6-7-3 face |
78 ! | U5, pressure onto the 3-7-8-4 face |
79 ! | U6, pressure onto the 4-8-5-1 face |
80 ! | U7, body force in the global X-direction |
81 ! | U8, body force in the global Y-direction |
82 ! | U9, body force in the global Z-direction |
83 ! | |
84 ! *--------------------------------------------------------------------*
85 
86  IMPLICIT DOUBLE PRECISION (a-h, o-z)
87 
88 ! Argument variables
89 
90  INTEGER :: numnp
91  INTEGER, DIMENSION(1:8,1:NumNp) :: nodes
92  REAL*8, DIMENSION(1:3*NumNp) :: rex
93  REAL*8, DIMENSION(1:3,1:NumNp) :: coords
94  INTEGER, DIMENSION(1:NumElVol) :: mapsfelvolel
95  dimension jdltyp(mdload,1:2),adlmag(mdload)
96 
97  INTEGER :: idvolel
98 
99 ! Local variables
100 
101  REAL*8, DIMENSION(1:3) :: veca, vecb, outwd_surface_normal, &
102  half_normal_veca, half_normal_vecb
103 
104  INTEGER :: iload
105  INTEGER :: ltype
106  INTEGER :: nd1, nd2, nd3, nd4, nd5, nd6, nd7, nd8
107 
108  REAL*8, DIMENSION(1:24) :: dlvec
109  INTEGER :: icnt
110 
111 ! Data statements
112 
113  DATA eigth /0.125000000000000/
114 
115 ! Compute external loads and add to the load vector
116 
117  DO iload = 1, ! LOOP OVER LOADS
118 
119  dlvec(:) = 0.d0
120 
121  idvolel = mapsfelvolel(i) ! Volumetric Element Id
122  ltype = jdltyp(iload,2) ! get load type
123 
124  nd1 = nodes(1,idvolel)
125  nd2 = nodes(2,idvolel)
126  nd3 = nodes(3,idvolel)
127  nd4 = nodes(4,idvolel)
128  nd5 = nodes(5,idvolel)
129  nd6 = nodes(6,idvolel)
130  nd7 = nodes(7,idvolel)
131  nd8 = nodes(8,idvolel)
132 
133  IF(ltype.LT.0) THEN ! unsupported load type (UnNU)
134  WRITE(6,1100)
135  stop
136  END IF
137 
138 ! branch depending on load type
139  SELECT CASE(ltype)
140 
141  CASE (100) ! pressure on Face 1
142 
143  veca(1:3) = coords(1:3,nd4) - coords(1:3,nd1)
144  vecb(1:3) = coords(1:3,nd2) - coords(1:3,nd1)
145 
146  half_normal_veca(1) = veca(2)*vecb(3) - vecb(2)*veca(3)
147  half_normal_veca(2) = veca(3)*vecb(1) - vecb(3)*veca(1)
148  half_normal_veca(3) = veca(1)*vecb(2) - vecb(1)*veca(2)
149 
150  veca(1:3) = coords(1:3,nd2) - coords(1:3,nd3)
151  vecb(1:3) = coords(1:3,nd4) - coords(1:3,nd3)
152 
153  half_normal_vecb(1) = veca(2)*vecb(3) - vecb(2)*veca(3)
154  half_normal_vecb(2) = veca(3)*vecb(1) - vecb(3)*veca(1)
155  half_normal_vecb(3) = veca(1)*vecb(2) - vecb(1)*veca(2)
156 
157  outwd_surface_normal = half_normal_veca + half_normal_vecb
158 
159  pmag = -eigth * adlmag(iload)
160 
161  dlvec(1:3) = pmag * outwd_surface_normal(1:3) ! node 1
162  dlvec(4:6) = pmag * outwd_surface_normal(1:3) ! node 2
163  dlvec(7:9) = pmag * outwd_surface_normal(1:3) ! node 3
164  dlvec(10:12) = pmag * outwd_surface_normal(1:3) ! node 4
165 
166  CASE(200) ! pressure on Face 2
167 
168  veca(1:3) = coords(1:3,nd6) - coords(1:3,nd5)
169  vecb(1:3) = coords(1:3,nd8) - coords(1:3,nd5)
170 
171  half_normal_veca(1) = veca(2)*vecb(3) - vecb(2)*veca(3)
172  half_normal_veca(2) = veca(3)*vecb(1) - vecb(3)*veca(1)
173  half_normal_veca(3) = veca(1)*vecb(2) - vecb(1)*veca(2)
174 
175  veca(1:3) = coords(1:3,nd8) - coords(1:3,nd7)
176  vecb(1:3) = coords(1:3,nd6) - coords(1:3,nd7)
177 
178  half_normal_vecb(1) = veca(2)*vecb(3) - vecb(2)*veca(3)
179  half_normal_vecb(2) = veca(3)*vecb(1) - vecb(3)*veca(1)
180  half_normal_vecb(3) = veca(1)*vecb(2) - vecb(1)*veca(2)
181 
182  outwd_surface_normal = half_normal_veca + half_normal_vecb
183 
184  pmag = -eigth * adlmag(iload)
185 
186  dlvec(13:15) = pmag * outwd_surface_normal(1:3) ! node 5
187  dlvec(16:18) = pmag * outwd_surface_normal(1:3) ! node 6
188  dlvec(19:21) = pmag * outwd_surface_normal(1:3) ! node 7
189  dlvec(22:24) = pmag * outwd_surface_normal(1:3) ! node 8
190 
191  CASE(300) ! pressure on Face 3
192 
193  veca(1:3) = coords(1:3,nd2) - coords(1:3,nd1)
194  vecb(1:3) = coords(1:3,nd5) - coords(1:3,nd1)
195 
196  half_normal_veca(1) = veca(2)*vecb(3) - vecb(2)*veca(3)
197  half_normal_veca(2) = veca(3)*vecb(1) - vecb(3)*veca(1)
198  half_normal_veca(3) = veca(1)*vecb(2) - vecb(1)*veca(2)
199 
200  veca(1:3) = coords(1:3,nd5) - coords(1:3,nd6)
201  vecb(1:3) = coords(1:3,nd2) - coords(1:3,nd6)
202 
203  half_normal_vecb(1) = veca(2)*vecb(3) - vecb(2)*veca(3)
204  half_normal_vecb(2) = veca(3)*vecb(1) - vecb(3)*veca(1)
205  half_normal_vecb(3) = veca(1)*vecb(2) - vecb(1)*veca(2)
206 
207  outwd_surface_normal = half_normal_veca + half_normal_vecb
208 
209  pmag = -eigth * adlmag(iload)
210 
211  dlvec(1:3) = pmag * outwd_surface_normal(1:3) ! node 1
212  dlvec(4:6) = pmag * outwd_surface_normal(1:3) ! node 2
213  dlvec(13:15) = pmag * outwd_surface_normal(1:3) ! node 5
214  dlvec(16:18) = pmag * outwd_surface_normal(1:3) ! node 6
215 
216  CASE(400) ! pressure on Face 4
217 
218  veca(1:3) = coords(1:3,nd3) - coords(1:3,nd2)
219  vecb(1:3) = coords(1:3,nd6) - coords(1:3,nd2)
220 
221  half_normal_veca(1) = veca(2)*vecb(3) - vecb(2)*veca(3)
222  half_normal_veca(2) = veca(3)*vecb(1) - vecb(3)*veca(1)
223  half_normal_veca(3) = veca(1)*vecb(2) - vecb(1)*veca(2)
224 
225  veca(1:3) = coords(1:3,nd6) - coords(1:3,nd7)
226  vecb(1:3) = coords(1:3,nd3) - coords(1:3,nd7)
227 
228  half_normal_vecb(1) = veca(2)*vecb(3) - vecb(2)*veca(3)
229  half_normal_vecb(2) = veca(3)*vecb(1) - vecb(3)*veca(1)
230  half_normal_vecb(3) = veca(1)*vecb(2) - vecb(1)*veca(2)
231 
232  outwd_surface_normal = half_normal_veca + half_normal_vecb
233  pmag = -eigth * adlmag(iload)
234 
235  dlvec(4:6) = pmag * outwd_surface_normal(1:3) ! node 2
236  dlvec(7:9) = pmag * outwd_surface_normal(1:3) ! node 3
237  dlvec(16:18) = pmag * outwd_surface_normal(1:3) ! node 6
238  dlvec(19:21) = pmag * outwd_surface_normal(1:3) ! node 7
239 
240 
241  CASE(500) ! pressure on Face 5
242 
243  veca(1:3) = coords(1:3,nd4) - coords(1:3,nd3)
244  vecb(1:3) = coords(1:3,nd7) - coords(1:3,nd3)
245 
246  half_normal_veca(1) = veca(2)*vecb(3) - vecb(2)*veca(3)
247  half_normal_veca(2) = veca(3)*vecb(1) - vecb(3)*veca(1)
248  half_normal_veca(3) = veca(1)*vecb(2) - vecb(1)*veca(2)
249 
250  veca(1:3) = coords(1:3,nd7) - coords(1:3,nd8)
251  vecb(1:3) = coords(1:3,nd4) - coords(1:3,nd8)
252 
253  half_normal_vecb(1) = veca(2)*vecb(3) - vecb(2)*veca(3)
254  half_normal_vecb(2) = veca(3)*vecb(1) - vecb(3)*veca(1)
255  half_normal_vecb(3) = veca(1)*vecb(2) - vecb(1)*veca(2)
256 
257  outwd_surface_normal = half_normal_veca + half_normal_vecb
258  pmag = -eigth * adlmag(iload)
259 
260  dlvec(7:9) = pmag * outwd_surface_normal(1:3) ! node 3
261  dlvec(10:12) = pmag * outwd_surface_normal(1:3) ! node 4
262  dlvec(19:21) = pmag * outwd_surface_normal(1:3) ! node 7
263  dlvec(22:24) = pmag * outwd_surface_normal(1:3) ! node 8
264 
265  CASE(600) ! pressure on Face 6
266 
267  veca(1:3) = coords(1:3,nd5) - coords(1:3,nd1)
268  vecb(1:3) = coords(1:3,nd4) - coords(1:3,nd1)
269 
270  half_normal_veca(1) = veca(2)*vecb(3) - vecb(2)*veca(3)
271  half_normal_veca(2) = veca(3)*vecb(1) - vecb(3)*veca(1)
272  half_normal_veca(3) = veca(1)*vecb(2) - vecb(1)*veca(2)
273 
274  veca(1:3) = coords(1:3,nd4) - coords(1:3,nd8)
275  vecb(1:3) = coords(1:3,nd5) - coords(1:3,nd8)
276 
277  half_normal_vecb(1) = veca(2)*vecb(3) - vecb(2)*veca(3)
278  half_normal_vecb(2) = veca(3)*vecb(1) - vecb(3)*veca(1)
279  half_normal_vecb(3) = veca(1)*vecb(2) - vecb(1)*veca(2)
280 
281  outwd_surface_normal = half_normal_veca + half_normal_vecb
282 
283  pmag = -eigth * adlmag(iload)
284 
285  dlvec(1:3) = pmag * outwd_surface_normal(1:3) ! node 1
286  dlvec(10:12) = pmag * outwd_surface_normal(1:3) ! node 4
287  dlvec(13:15) = pmag * outwd_surface_normal(1:3) ! node 5
288  dlvec(22:24) = pmag * outwd_surface_normal(1:3) ! node 8
289 
290  CASE(700) ! body force in X-direction
291 
292  CASE(800) ! body force in Y-direction
293 
294  CASE(900)! body force in Z-direction
295 
296 
297 
298  END SELECT
299 
300  icnt = 1
301  DO i = 1, 8
302  nd = nodes(i,idvolel)
303  rex(3*nd-2) = rex(3*nd-2) + dlvec(icnt)
304  rex(3*nd-1) = rex(3*nd-1) + dlvec(icnt+1)
305  rex(3*nd ) = rex(3*nd ) + dlvec(icnt+2)
306  icnt = icnt + 3
307  ENDDO
308 
309  END DO ! END LOOP OVER LOADS
310 
311 1100 FORMAT (/,2x,'>>> Non-uniform surface tractions are not', &
312  /,2x,' supported (i.e. UnNU). Aborting !')
313 
314 END SUBROUTINE get_external_loads
315 
get_external_loads
subroutine get_external_loads(coords, jdltyp, adlmag, mdload, ndload, Rex, nodes, NumNp, MapSElVolEl, NumElVol)
Definition: get_external_loads.f90:53
Connectivity::dimension
int dimension() const
Get the dimension of the mesh.
Definition: Connectivity.h:125
z
void int int int REAL REAL REAL * z
Definition: write.cpp:76
i
blockLoc i
Definition: read.cpp:79
x
void int int REAL * x
Definition: read.cpp:74
Line_2::a
RT a() const
Definition: Line_2.h:140