/[escript]/trunk/doc/cookbook/example09.tex
ViewVC logotype

Diff of /trunk/doc/cookbook/example09.tex

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3152 by jfenwick, Mon Jul 26 06:10:50 2010 UTC revision 3153 by ahallam, Sun Sep 5 01:31:49 2010 UTC
# Line 12  Line 12 
12  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
13    
14  \section{3D pycad}  \section{3D pycad}
15  \sslist{example09n.py}  \sslist{example09m.py}
16    This example explains how to build a 3D layered domain using pycad. As
17  This example explains how to build a 3D layered model using pycad. There are a  simple as this example sounds, there are a few import concepts that one must
18  few import concepts that one must remember when using pycad to build a layered  remember  so that the model will function correctly.
 model that will function correctly.  
19  \begin{itemize}  \begin{itemize}
20    \item There must be no duplication of any geometric features whether they be    \item There must be no duplication of any geometric features whether they are
21    points, lines, loops, surfaces or volumes.    points, lines, loops, surfaces or volumes.
22    \item All objects with dimensions greater then a line have a normal defined by    \item All objects with dimensions greater then a line have a normal defined by
23    the right hand rull (RHR). It is important to consider which direction a    the right hand rull (RHR). It is important to consider which direction a
# Line 38  import os Line 37  import os
37  \end{python}  \end{python}
38  After carrying out some routine checks and setting the \verb!save_path! we then  After carrying out some routine checks and setting the \verb!save_path! we then
39  specify the parameters of the model. This model will be 2000 by 2000 meters on  specify the parameters of the model. This model will be 2000 by 2000 meters on
40  the surface and extend to a depth of 500 meters. An interface of boundary  the surface and extend to a depth of 500 meters. An interface or boundary
41  between two layers will be created at half the total depth or 250 meters. This  between two layers will be created at half the total depth or 250 meters. This
42  type of model is known as a horizontally layered model or a layer cake model.  type of model is known as a horizontally layered model or a layer cake model.
43  \begin{python}  \begin{python}
# Line 84  linhe_ar.append(Line(ip3,ip0)) Line 83  linhe_ar.append(Line(ip3,ip0))
83  Consider now the sides of the domain. One could specify the whole side using the  Consider now the sides of the domain. One could specify the whole side using the
84  points first defined for the top and bottom layer. This would specify the whole  points first defined for the top and bottom layer. This would specify the whole
85  domain as one volume. However, there is an interface and we wish to define each  domain as one volume. However, there is an interface and we wish to define each
86  layer individually. Therefore there will be 8 surfaces on the sides of our  layer individually. Therefore, there will be 8 surfaces on the sides of our
87  domain. We can do this operation quite simply using the points and lines that we  domain. We can do this operation quite simply using the points and lines that we
88  had defined previously. First loops are created and then surfaces making sure to  had defined previously. First loops are created and then surfaces making sure to
89  keep a normal for each layer which is consistent with upper and lower surfaces  keep a normal for each layer which is consistent with upper and lower surfaces
90  of the layer. For example all normals must face outwards from or inwards towards  of the layer. For example, all surface normals must face outwards from or
91   the centre of the volume.  inwards towards the centre of the volume.
92  \begin{python}  \begin{python}
93  cintfa_ar=[]; cintfb_ar=[] #curveloops for above and below interface on sides  cintfa_ar=[]; cintfb_ar=[] #curveloops for above and below interface on sides
94  cintfa_ar.append(CurveLoop(linte_ar[0],linhe_ar[0],-linte_ar[2],-l01))  cintfa_ar.append(CurveLoop(linte_ar[0],linhe_ar[0],-linte_ar[2],-l01))
# Line 102  sintfb_ar=[PlaneSurface(cintfb_ar[i]) fo Line 101  sintfb_ar=[PlaneSurface(cintfb_ar[i]) fo
101  sintf=PlaneSurface(CurveLoop(*tuple(linhe_ar)))  sintf=PlaneSurface(CurveLoop(*tuple(linhe_ar)))
102  \end{python}  \end{python}
103  Assuming all is well with the normals, the volumes can be created from our  Assuming all is well with the normals, the volumes can be created from our
104  surface arrays. Note the use here of the \verb!*tuple*! function. This allows us  surface arrays. Note the use here of the \verb!*tuple! function. This allows us
105  to pass an list array as an argument list to a function. It must be placed at  to pass an list array as an argument list to a function. It must be placed at
106  the end of the function arguments and there cannot be more than one per function  the end of the function arguments and there cannot be more than one per function
107  call.  call.
# Line 165  from pycad. Line 164  from pycad.
164  from esys.pycad import layer_cake  from esys.pycad import layer_cake
165  intfaces=[10,30,50,55,80,100,200,250,400,500]  intfaces=[10,30,50,55,80,100,200,250,400,500]
166    
167  cmplx_domain=layer_cake.LayerCake(xwidth,ywidth,intfaces,200.0)  domaindes=Design(dim=3,element_size=element_size,order=2)
168    cmplx_domain=layer_cake.LayerCake(domaindes,xwidth,ywidth,intfaces)
169  cmplx_domain.setScriptFileName(os.path.join(save_path,"example09lc.geo"))  cmplx_domain.setScriptFileName(os.path.join(save_path,"example09lc.geo"))
170  cmplx_domain.setMeshFileName(os.path.join(save_path,"example09lc.msh"))  cmplx_domain.setMeshFileName(os.path.join(save_path,"example09lc.msh"))
171  dcmplx=MakeDomain(cmplx_domain)  dcmplx=MakeDomain(cmplx_domain)
# Line 187  directionality of your primatives when c Line 187  directionality of your primatives when c
187  one cannot establist the tangent of a line or curveloop, or the normal of a  one cannot establist the tangent of a line or curveloop, or the normal of a
188  surface. These values can be checked by importing the geometry to gmesh and  surface. These values can be checked by importing the geometry to gmesh and
189  applying the appropriate options.  applying the appropriate options.
190    
191    \section{3D Seismic Wave Propagation}
192    Adding an extra dimension to our wave equation solution script should be
193    relatively simple. Apart from the changes to the domain and therefore the
194    parameters of the model, there are only a few minor things which must be
195    modified to make the solution appropriate for 3d modelling.
196    
197    \section{Applying a function to a domain tag}
198    \sslist{example09b.py}
199    To apply a function or data object to a domain requires a two step process. The
200    first step is to create a data object with a on/off mask based upon the tagged
201    value. This is quite simple and can be achieved by using a scalar data object
202    based upon the domain. In this case we are using the \verb!FunctionOnBoundary!
203    function space because the tagged value \verb!'stop'! is effectively a specific
204    subsurface of the boundary of the whole domain.
205    \begin{python}
206    stop=Scalar(0.0,FunctionOnBoundary(domain))
207    stop.setTaggedValue("stop",1.0)
208    \end{python}
209    Now the data object \verb|stop| has the value 1.0 on the surface
210    \verb!'stop'! and zero everywhere else.
211    %
212     To apply our function to the boundary only on \verb|stop| we now
213     mulitply it by \verb|stop|
214    \begin{python}
215    xb=FunctionOnBoundary(domain).getX()
216    yx=(cos(length(xb-xc)*3.1415/src_length)+1)*whereNegative(length(xb-xc)-src_length)
217    src_dir=numpy.array([0.,0.,1.0]) # defines direction of point source as down
218    mypde.setValue(y=source[0]*yx*src_dir*stop) #set the source as a function on the boundary
219    \end{python}
220    
221    \section{Mayavi2 with 3D data.}
222    There are a variety of visualisation options available when using VTK data. The
223    types of visualisation are often data and interpretation specific and thus
224    consideration must be given to the type of output saved to file. Whether that is
225    scalar, vector or tensor data.
226    
227    \begin{figure}[htp]
228    \centering
229    \begin{subfigure}[Example 9 at time step 201.]
230        {\label{fig:ex9b 201}
231        \includegraphics[width=0.45\textwidth]{figures/ex09b00201.png}}
232     \end{subfigure}
233     \begin{subfigure}[Example 9 at time step 341.]
234        {\label{fig:ex9b 201}
235        \includegraphics[width=0.45\textwidth]{figures/ex09b00341.png}}
236     \end{subfigure}\\
237     \begin{subfigure}[Example 9 at time step 440.]
238        {\label{fig:ex9b 201}
239        \includegraphics[width=0.45\textwidth]{figures/ex09b00440.png}}
240     \end{subfigure}
241    \end{figure}

Legend:
Removed from v.3152  
changed lines
  Added in v.3153

  ViewVC Help
Powered by ViewVC 1.1.26