[cig-commits] r19964 - in seismo/3D/SPECFEM3D_GEOTECH/trunk: . doc src
homnath at geodynamics.org
homnath at geodynamics.org
Mon Apr 23 09:23:34 PDT 2012
Author: homnath
Date: 2012-04-23 09:23:33 -0700 (Mon, 23 Apr 2012)
New Revision: 19964
Modified:
seismo/3D/SPECFEM3D_GEOTECH/trunk/TODO
seismo/3D/SPECFEM3D_GEOTECH/trunk/doc/cmake.pdf
seismo/3D/SPECFEM3D_GEOTECH/trunk/doc/cover.pdf
seismo/3D/SPECFEM3D_GEOTECH/trunk/doc/manual_SPECFEM3D_GEOTECH.bib
seismo/3D/SPECFEM3D_GEOTECH/trunk/doc/manual_SPECFEM3D_GEOTECH.tex
seismo/3D/SPECFEM3D_GEOTECH/trunk/src/Makefile
seismo/3D/SPECFEM3D_GEOTECH/trunk/src/apply_traction.f90
seismo/3D/SPECFEM3D_GEOTECH/trunk/src/compute_bmat.f90
seismo/3D/SPECFEM3D_GEOTECH/trunk/src/compute_cmat.f90
seismo/3D/SPECFEM3D_GEOTECH/trunk/src/excavation.f90
seismo/3D/SPECFEM3D_GEOTECH/trunk/src/gll_library.f90
seismo/3D/SPECFEM3D_GEOTECH/trunk/src/global.f90
seismo/3D/SPECFEM3D_GEOTECH/trunk/src/math_library.f90
seismo/3D/SPECFEM3D_GEOTECH/trunk/src/plastic_library.f90
seismo/3D/SPECFEM3D_GEOTECH/trunk/src/preprocess.f90
seismo/3D/SPECFEM3D_GEOTECH/trunk/src/read_input.f90
seismo/3D/SPECFEM3D_GEOTECH/trunk/src/semexcav3d.F90
seismo/3D/SPECFEM3D_GEOTECH/trunk/src/semgeotech.F90
seismo/3D/SPECFEM3D_GEOTECH/trunk/src/semslope3d.F90
seismo/3D/SPECFEM3D_GEOTECH/trunk/src/string_library.f90
Log:
Minor update
Modified: seismo/3D/SPECFEM3D_GEOTECH/trunk/TODO
===================================================================
--- seismo/3D/SPECFEM3D_GEOTECH/trunk/TODO 2012-04-21 00:26:41 UTC (rev 19963)
+++ seismo/3D/SPECFEM3D_GEOTECH/trunk/TODO 2012-04-23 16:23:33 UTC (rev 19964)
@@ -1,6 +1,7 @@
REVISION:
HNG, Jul 08,2011
TODO LIST
+ - add unrecognized variables feature in input
- automatic factor of safey determination
- need to consider ghost partitions for the averaging nodal quantities
- block visualization so that interior part can be clearly visualized
Modified: seismo/3D/SPECFEM3D_GEOTECH/trunk/doc/cmake.pdf
===================================================================
(Binary files differ)
Modified: seismo/3D/SPECFEM3D_GEOTECH/trunk/doc/cover.pdf
===================================================================
(Binary files differ)
Modified: seismo/3D/SPECFEM3D_GEOTECH/trunk/doc/manual_SPECFEM3D_GEOTECH.bib
===================================================================
--- seismo/3D/SPECFEM3D_GEOTECH/trunk/doc/manual_SPECFEM3D_GEOTECH.bib 2012-04-21 00:26:41 UTC (rev 19963)
+++ seismo/3D/SPECFEM3D_GEOTECH/trunk/doc/manual_SPECFEM3D_GEOTECH.bib 2012-04-23 16:23:33 UTC (rev 19964)
@@ -79,7 +79,7 @@
owner = {homnath},
timestamp = {2011.05.27},
url = {http://cubit.sandia.gov/}
-}
+}
@MANUAL{ensight2008,
@@ -92,7 +92,7 @@
note = {[Online; accessed 11-July-2011]},
owner = {homnath},
timestamp = {2011.07.11}
-}
+}
@ARTICLE{faccioli1997,
author = {Faccioli, E. and Maggio, F. and Paolucci, R. and Quarteroni, A.},
@@ -106,7 +106,7 @@
issue = {3},
keyword = {Earth and Environmental Science},
publisher = {Springer Netherlands}
-}
+}
@ARTICLE{geuzaine2009,
author = {C. Geuzaine and J. F. Remacle},
@@ -121,28 +121,29 @@
timestamp = {2011.01.23}
}
- at ARTICLE{gharti2011,
+ at ARTICLE{gharti2012b,
author = {Gharti, H. N. and Komatitsch, D. and Oye, V. and Martin, R. and Tromp,
J.},
title = {Application of an elastoplastic spectral-element method to {3D} slope
stability analysis},
journal = ijnme,
- year = {2011},
- volume = {submitted},
+ year = {2012},
+ volume = {in press},
+ doi = {http://dx.doi.org/10.1002/nme.3374.},
owner = {homnath},
timestamp = {2011.06.01}
}
- at ARTICLE{gharti2011b,
- author = {Gharti, H. N. and Oye, V. and Komatitsch, D. and Tromp, J.},
- title = {Simulation of multistage excavation based on a {3D} spectral-element
+ at ARTICLE{gharti2012a,
+ author = {Gharti, H. N. and Oye, V. and Komatitsch, D. and Tromp, J.},
+ title = {Simulation of multistage excavation based on a 3D spectral-element
method},
- journal = {Computers \& Structures},
- year = {2011},
- volume = {submitted},
- owner = {homnath},
- timestamp = {2011.06.01}
-}
+ journal = cs,
+ year = {2012},
+ volume = {100--101},
+ pages = {54--69},
+ doi = {10.1016/j.compstruc.2012.03.005}
+}
@BOOK{gropp1994,
title = {Using {MPI}, portable parallel programming with the {M}essage-{P}assing
@@ -164,7 +165,7 @@
number = {2},
doi = {10.1016/0045-7825(83)90115-9},
issn = {0045-7825},
-}
+}
@ARTICLE{khan1996,
author = {Khan, A. I. and Topping, B. H. V.},
@@ -188,7 +189,7 @@
number = {1},
doi = {10.1016/0045-7825(87)90182-4},
issn = {0045-7825},
-}
+}
@ARTICLE{komatitsch1998,
author = {Komatitsch, D. and Vilotte, J. P.},
@@ -199,7 +200,7 @@
volume = {88},
pages = {368--392},
number = {2},
-}
+}
@ARTICLE{komatitsch1999,
author = {Komatitsch, D. and Tromp, J.},
@@ -214,7 +215,7 @@
modelling, seismic wave propagation, topography.},
owner = {homnath},
timestamp = {2009.07.22}
-}
+}
@TECHREPORT{larsen1995,
author = {Larsen, S. and Schultz, C. A.},
@@ -223,7 +224,7 @@
year = {1995},
owner = {homnath},
timestamp = {2008.04.08}
-}
+}
@ARTICLE{law1986,
author = {Law, K. H.},
@@ -235,7 +236,7 @@
number = {6},
doi = {10.1016/0045-7949(86)90254-3},
issn = {0045-7949},
-}
+}
@BOOK{pacheco1997,
title = {Parallel Programming with {MPI}},
@@ -260,7 +261,7 @@
keywords = {Tomography, Interferometry, Computational seismology, Wave propagation},
publisher = {Blackwell Publishing Ltd},
url = {http://dx.doi.org/10.1111/j.1365-246X.2011.05044.x}
-}
+}
@ARTICLE{patera1984,
author = {Patera, A. T.},
Modified: seismo/3D/SPECFEM3D_GEOTECH/trunk/doc/manual_SPECFEM3D_GEOTECH.tex
===================================================================
--- seismo/3D/SPECFEM3D_GEOTECH/trunk/doc/manual_SPECFEM3D_GEOTECH.tex 2012-04-21 00:26:41 UTC (rev 19963)
+++ seismo/3D/SPECFEM3D_GEOTECH/trunk/doc/manual_SPECFEM3D_GEOTECH.tex 2012-04-23 16:23:33 UTC (rev 19964)
@@ -51,19 +51,12 @@
% \includegraphics[width=\paperwidth]{cover}
%\end{textblock*}
-%\begin{figure}[H]
-%\centering
-%\noindent\includegraphics[width=0.8\paperwidth]{cover}
-%\noindent\includegraphics[scale=1.0]{cover}
-%\end{figure}
-
\thispagestyle{empty} % no page number
\title{\textbf{\packver \\
User Manual}}
-
-\author{Hom Nath Gharti, NORSAR, Norway \\
-Dimitri Komatitsch, University of Toulouse, France \\
+\author{Hom Nath Gharti\footnote{Previously at: NORSAR, Norway}, Princeton University, USA \\
+Dimitri Komatitsch, Aix-Marseille University, France\\
Volker Oye, NORSAR, Norway \\
Roland Martin, University of Toulouse, France \\
Jeroen Tromp, Princeton University, USA}
@@ -108,13 +101,13 @@
\chapter{Introduction}
\section{Background}
-\pack\ is a free and open-source command-driven software for 3D slope stability analysis~\citep[For more detail see][]{gharti2011} and simulation of 3D multistage excavation~\citep[For more detail see][]{gharti2011b} based on the spectral-element method~\citep[e.g.,][]{patera1984,canuto1988,seriani1994,faccioli1997,komatitsch1998,komatitsch1999,peter2011}. The software can run on a single processor as well as multi-core machines or large clusters. It is written mainly in FORTRAN 90, and parallelized using
-MPI~\citep{gropp1994,pacheco1997} based on domain decomposition. For the domain decomposition, an open-source graph partitioning library SCOTCH~\citep{pellegrini1996} is used. The element-by-element preconditioned conjugate-gradient method~\citep[e.g.,][]{hughes1983,law1986,king1987,barragy1988} is implemented to solve the linear equations. For elastoplastic failure,
-Mohr-coulomb failure criterion is used with viscoplastic strain method~ \citep{zienkiewicz1974}.\\
+\pack\ is a free and open-source command-driven software for 3D slope stability analysis~\citep[for more details see][]{gharti2012b} and simulation of 3D multistage excavation~\citep[for more details see][]{gharti2012a} based on the spectral-element method~\citep[e.g.,][]{patera1984,canuto1988,seriani1994,faccioli1997,komatitsch1998,komatitsch1999,peter2011}. The slope stability and the excavation routines were originally started from the routines found in the book ``Programming the finite element method''~\citep{smith2004}. The software can run on a single processor as well as multi-core machines or large clusters. It is written mainly in FORTRAN 90, and parallelized using
+MPI~\citep{gropp1994,pacheco1997} based on domain decomposition. For the domain decomposition, the open-source graph partitioning library SCOTCH~\citep{pellegrini1996} is used. The element-by-element preconditioned conjugate-gradient method~\citep[e.g.,][]{hughes1983,law1986,king1987,barragy1988} is implemented to solve the linear equations. For elastoplastic failure,
+a Mohr-coulomb failure criterion is used with a viscoplastic strain method~ \citep{zienkiewicz1974}.\\
-This program does not automatically determine the factor of safety of the slope stability. Simulation can be performed for a series of safety factors. After plotting the safety factor vs maximum displacement curve, one can determine the factor of safety of the given slope. Although, the software is optimized for slope stability analysis and multistage excavation, other relevant simulations of static problems in solid (geo)mechanics can also be performed with this software.\\
+This program does not automatically determine the factor of safety of slope stability. Simulations can be performed for a series of safety factors. After plotting the safety factor verses maximum displacement curve, one can determine the factor of safety of the given slope. Although the software is optimized for slope stability analysis and multistage excavation, other relevant simulations of quasistatic problems in solid (geo)mechanics can also be performed with this software.\\
-The software currently does not include the inbuilt mesher. Existing tools such as Gmsh~\citep{geuzaine2009}, CUBIT~\citep{cubit2011}, TrueGrid~\citep{truegrid2006}, etc. can be used for the hexahedral meshing, and the resulting mesh file can be converted to the input files required by the \pack. Output data can be visualized and processed using an open-source visualization application ParaView (\texttt{www.paraview.org}).
+The software currently does not include an inbuilt mesher. Existing tools, such as Gmsh~\citep{geuzaine2009}, CUBIT~\citep{cubit2011}, TrueGrid~\citep{truegrid2006}, etc., can be used for hexahedral meshing, and the resulting mesh file can be converted to the input files required by \pack. Output data can be visualized and processed using the open-source visualization application ParaView (\texttt{www.paraview.org}).
\section{Status summary}
\begin{desclist}{Pseudo-static earthquake loading}
@@ -127,23 +120,21 @@
\item[Automatic factor of safety] : No
\end{desclist}
-\section*{Revision}
+\section*{Revisions}
-HNG, Sep 08, 2011; HNG, Jul 12, 2011; HNG, May 20, 2011; HNG, Jan 17, 2011
+HNG, Jan 12, 2012; HNG, Sep 08, 2011; HNG, Jul 12, 2011; HNG, May 20, 2011; HNG, Jan 17, 2011
\chapter{Getting started}
\section{Package structure}
-Original \pack\ package comes in a single compressed file \linebreak\texttt{\pack.tar.gz}, which can be extracted using \texttt{tar} command:\\
+The original \pack\ package comes in a single compressed file \linebreak\texttt{\pack.tar.gz}, which can be extracted using \texttt{tar} command:\\
\texttt{tar -zxvf \pack.tar.gz}\\
-Or\\
+or using, for example, \texttt{7-zip (www.7-zip.org)} under WINDOWS. The package has the following structure:\\
-using, for example, \texttt{7-zip (www.7-zip.org)} in WINDOWS. The package has a following structure.\\
-
\texttt{\pack/}
\begin{adescription}{~~CMakeLists.txt}
\item[~~COPYING] : License.
@@ -159,17 +150,17 @@
\section{Prerequisites}
\begin{itemize}[-]
- \item \underline{CMake build system}. The CMake version >= 2.8.4 is necessary to configure the software. It is free and open-source, and can be downloaded from \texttt{www.cmake.org}.
- \item \underline{Make utility}. The make utility is necessary to build the software using Makefile. This utility is usually installed by default in most of the LINUX systems. In WINDOWS, one can use Cygwin (\texttt{www.cygwin.com}) or MinGW (\texttt{www.mingw.org}) to install the make utility.
- \item \underline{A recent FORTRAN compiler}. The software is written mainly in FORTRAN 90, but it also uses a few FORTRAN 2003 features (e.g., streaming IO). These features are already available in most of the FORTRAN compilers, e.g., gfortran version >= 4.2 (\texttt{gcc.gnu.org/wiki/GFortran}) and g95 (\texttt{www.g95.org}).
+ \item \underline{CMake build system}. The CMake version $\ge$ 2.8.4 is necessary to configure the software. It is free and open-source, and can be downloaded from \texttt{www.cmake.org}.
+ \item \underline{Make utility}. The make utility is necessary to build the software using Makefile. This utility is usually installed by default in most LINUX systems. Under WINDOWS, one can use Cygwin (\texttt{www.cygwin.com}) or MinGW (\texttt{www.mingw.org}) to install the make utility.
+ \item \underline{A recent FORTRAN compiler}. The software is written mainly in FORTRAN 90, but it also uses a few FORTRAN 2003 features (e.g., streaming IO). These features are already available in most of the FORTRAN compilers, e.g., gfortran version $\ge$ 4.2 (\texttt{gcc.gnu.org/wiki/GFortran}) and g95 (\texttt{www.g95.org}).
\end{itemize}
Following libraries are necessary for parallel processing.
\begin{itemize}[-]
- \item \underline{A recent MPI library}. It should be built with same FORTRAN compiler which will be used to compile the software. Please see \texttt{www.open-mpi.org} or \linebreak\texttt{www.mcs.anl.gov/research/projects/mpich2} for detail on how to install MPI library and how to run MPI programs.
- \item \underline{SCOTCH graph partitioning library}. This library should be compiled with same\linebreak FORTRAN compiler which will be used to compile the software. Please see\linebreak \texttt{www.labri.fr/perso/pelegrin/scotch} for detail on how to install SCOTCH. Version 5.1.7 was successfully tested with \pack.
+ \item \underline{A recent MPI library}. It should be built with the same FORTRAN compiler used to compile the software. Please see \texttt{www.open-mpi.org} or \linebreak\texttt{www.mcs.anl.gov/research/projects/mpich2} for details on how to install MPI library and how to run MPI programs.
+ \item \underline{SCOTCH graph partitioning library}. This library should be compiled with the same\linebreak FORTRAN compiler used to compile the software. Please see\linebreak \texttt{www.labri.fr/perso/pelegrin/scotch} for details on how to install SCOTCH.
\end{itemize}
- Finally, following compiler is necessary to build the documentation (this file).
+ Finally, the following compiler is necessary to build the documentation (this file):
\begin{itemize}[-]
\item \underline{\LaTeX\ compiler}. This is necessary to compile the documentation files.
\end{itemize}
@@ -177,16 +168,16 @@
\section{Configure}
\label{sec:configure}
-Software package \pack\ is configured using CMake, and the package uses out-of-source build. Hence, \underline{DO NOT} build in the same source directory. Let's say the full path to the package (source directory) is \texttt{\$HOME/download/\pack}.
+Software package \pack\ is configured using CMake, and the package uses an out-of-source build. Hence, \underline{DO NOT} build in the same source directory. Let's say the full path to the package (source directory) is \texttt{\$HOME/download/\pack}.
\begin{itemize}
\item Create a separate build directory, e.g.,\\
\texttt{mkdir \$HOME/work/\pack}
-\item Go to build directory \\
+\item Go to the build directory \\
\texttt{cd \$HOME/work/\pack}
-\item Type cmake command \\
+\item Type the cmake command \\
\texttt{ccmake \$HOME/projects/\pack}
\end{itemize}
@@ -204,8 +195,8 @@
\item Configure (c key or Configure button)
\end{itemize}
-If WARNINGS or ERRORS occur, press e key (or OK button) to return to configuration. These steps have to be repeated until successful configuration. Then, press g key (or Generate button) to generate build files. Check carefully that all necessary variables are set properly. Unless configuration is successful, generate is not enabled. Sometimes, c key (or Configure button) has to be pressed repeatedly until generate is enabled. Initially, all variables may not be visible. To see all variables, toggle advanced mode pressing t key (or Advanced button).
-To set or change a variable, move the cursor to the variable and press Enter key. If the variable is a boolean (ON/OFF), it will flip the value on pressing the Enter key. If the variable is a string or a file, it can be edited. For more detail, please see the CMake documentation (\texttt{www.cmake.org}).
+If WARNINGS or ERRORS occur, press the e key (or the OK button) to return to configuration. These steps have to be repeated until successful configuration. Then, press the g key (or the Generate button) to generate build files. Check carefully that all necessary variables are set properly. Unless configuration is successful, generate is not enabled. Sometimes, the c key (or Configure button) has to be pressed repeatedly until generate is enabled. Initially, all variables may not be visible. To see all variables, toggle advanced mode by pressing the t key (or the Advanced button).
+To set or change a variable, move the cursor to the variable and press Enter key. If the variable is a boolean (ON/OFF), it will flip the value on pressing the Enter key. If the variable is a string or a file, it can be edited. For more details, please see the CMake documentation (\texttt{www.cmake.org}).
\\
Following are the main CMake variables for the \pack\ (See Figure~\ref{fig:cmake})
@@ -213,25 +204,25 @@
\colorbox{gray}{
\parbox{15.5cm}{
\begin{adescription}{BUILD\_UTILITIES\_EXODUS2SEM}
-\item[BUILD\_DOCUMENTATION] : If \texttt{ON}, user manual (this file) is created. The default is \texttt{OFF}.
-\item[BUILD\_PARTMESH] : If \texttt{ON}, \texttt{partmesh} program is built. The default is \texttt{OFF}. The \texttt{partmesh} program is necessary to partition the mesh for parallel processing.
-\item[BUILD\_UTILITIES\_EXODUS2SEM] : If \texttt{ON}, \texttt{exodus2sem} program is built. The default is \texttt{OFF}. The \texttt{exodus2sem} program convert exodus mesh file to input files required by the \pack\ package (see also Chapter~\ref{chap:utilities}).
-\item[BUILD\_UTILITIES\_WRITE\_SOS] : If \texttt{ON}, \texttt{write\_sos} program is built. The default is \texttt{OFF}. The \texttt{write\_sos} program writes a EnSight SOS file necessary for the parallel visualization (see also Chapter~\ref{chap:utilities}).
-\item[ENABLE\_MPI] : If \texttt{ON}, main parallel program \texttt{psemgeotech} is built otherwise main serial program \texttt{semgeotech} is built. The default is \texttt{OFF}.
+\item[BUILD\_DOCUMENTATION] : If \texttt{ON}, the user manual (this file) is created. The default is \texttt{OFF}.
+\item[BUILD\_PARTMESH] : If \texttt{ON}, the \texttt{partmesh} program is built. The default is \texttt{OFF}. The \texttt{partmesh} program is necessary to partition the mesh for parallel processing.
+\item[BUILD\_UTILITIES\_EXODUS2SEM] : If \texttt{ON}, the \texttt{exodus2sem} program is built. The default is \texttt{OFF}. The \texttt{exodus2sem} program convert exodus mesh file to input files required by the \pack\ package (see also Chapter~\ref{chap:utilities}).
+\item[BUILD\_UTILITIES\_WRITE\_SOS] : If \texttt{ON}, the \texttt{write\_sos} program is built. The default is \texttt{OFF}. The \texttt{write\_sos} program writes a EnSight SOS file necessary for the parallel visualization (see also Chapter~\ref{chap:utilities}).
+\item[ENABLE\_MPI] : If \texttt{ON}, the main parallel program \texttt{psemgeotech} is built otherwise main serial program \texttt{semgeotech} is built. The default is \texttt{OFF}.
\item[SCOTCH\_LIBRARY\_PATH] : This is required if \texttt{BUILD\_PARTMESH} is \texttt{ON}. If not found automatically, it can be set manually.
-\item[CMAKE\_Fortran\_COMPILER] : This defines the Fortran compiler. If not found automatically or automatically found compiler is not correct, it can be set manually.
+\item[CMAKE\_Fortran\_COMPILER] : This defines the Fortran compiler. If not found automatically or the automatically found compiler is not correct, it can be set manually.
\end{adescription}
}}\\
{\emph{Note 1: If}} \texttt{CMAKE\_Fortran\_COMPILER} {\emph{has to be changed, first change this and configure, and then change other variables if necessary and configure.}}\\
-{\emph{Note 2: Even if some of the above variables are set }} \texttt{ON}{\emph{, if appropriate working compilers are not found, corresponding variables are internally set}} \texttt{OFF} {\emph{with WARNING message.}}
+{\emph{Note 2: Even if some of the above variables are set }} \texttt{ON}{\emph{, if appropriate working compilers are not found, corresponding variables are internally set}} \texttt{OFF} {\emph{with WARNING messages.}}
\section{Compile}
\begin{itemize}[]
- \item Once configuration and generation are successful, necessary build files are created. Now to build the main program, type: \\
+ \item Once configuration and generation are successful, the necessary build files are created. Now to build the main program, type: \\
\texttt{make}
- \item On the multi-processor system (let's say eight processor), type:\\
+ \item On multi-processor systems (let's say eight processors), type:\\
\texttt{make -j 8}
\item To clean, type\\
@@ -242,7 +233,7 @@
\section{Run}
\subsubsection{Serial run}
\begin{itemize}[-]
-\item To run serial program, type \\
+\item To run the serial program, type \\
\texttt{./bin/semgeotech} \emph{input\_file\_name}
Example:
@@ -260,7 +251,7 @@
\texttt{./bin/partmesh ./input/validation1.psem}
-\item To run parallel program, type \\
+\item To run the parallel program, type \\
\texttt{mpirun -n} \emph{number\_of\_nodes} \texttt{./bin/psemgeotech} \emph{input\_file\_name} \\
OR
@@ -272,13 +263,13 @@
\texttt{mpirun -n 8 ./bin/psemgeotech ./input/validation1.psem}
\end{itemize}
-{\emph{Note: see Chapter~\ref{chap:input} for detail on input and input files. Try to run one or more examples included in}} \texttt{input/}{\emph{. By default, example files included in the package are not copied to build directory during build process. If necessary, copy files within}} \texttt{input/} {\emph{folder of source directory to the}} \texttt{input/} {\emph{folder of build directory.}}
+{\emph{Note: see Chapter~\ref{chap:input} for details on input and input files. Try to run one or more examples included in}} \texttt{input/}{\emph{. By default, example files included in the package are not copied to build directory during build process. If necessary, copy files within}} \texttt{input/} {\emph{folder of source directory to the}} \texttt{input/} {\emph{folder of build directory.}}
\chapter{Input}
\label{chap:input}
\section{Main input file}
-Main input file structure is motivated by the ``E3D''~\citep{larsen1995} software package. The main input file consists of legitimate input lines defined in the specified formats. Any number of blank lines or comment lines can be placed for user friendly input structure. The blank lines contain no or only white-space characters, and the comment lines contain "\#" as the first character. \\
+The main input file structure is motivated by the ``E3D''~\citep{larsen1995} software package. The main input file consists of legitimate input lines defined in the specified formats. Any number of blank lines or comment lines can be placed for user friendly input structure. The blank lines contain no or only white-space characters, and the comment lines contain "\#" as the first character. \\
Each legitimate input line consists of a line type, and list of arguments and corresponding values. All argument-value pair are separated by comma (,). If necessary, any legitimate input line can be continued to next line using FORTRAN 90 continuation character "\&" as an absolute last character of a line to be continued. Repetition of same line type is not allowed.\\
@@ -289,7 +280,7 @@
\texttt{preinfo: nproc=8, ngllx=3, nglly=3, ngllz=3, nenod=8, ngnod=8, \& \\
inp\_path=\sq{../input}, part\_path=\sq{../partition}, out\_path=\sq{../output/}}\\
-All legitimate input lines should be written in lower case. Line type and argument-value pairs must be separated by space. Each argument-value pair must be separated by comma(,) and space/s. No space/s are recommended before line type and in between argument name and "=" or "=" and argument value. If argument value is a string, the FORTRAN 90 string (i.e., enclosed within the single quotes) should be used, for example, \texttt{inp\_path=\sq{../input}}. If the argument value is a vector (i.e., multi-valued), a list of values separated by space (no comma!) shoud be used, e.g, \texttt{srf=1.0 1.2 1.3 1.4}.
+All legitimate input lines should be written in lower case. Line type and argument-value pairs must be separated by a space. Each argument-value pair must be separated by a comma(,) and a space/s. No space/s are recommended before line type and in between argument name and "=" or "=" and argument value. If argument value is a string, the FORTRAN 90 string (i.e., enclosed within single quotes) should be used, for example, \texttt{inp\_path=\sq{../input}}. If the argument value is a vector (i.e., multi-valued), a list of values separated by space (no comma!) shoud be used, e.g, \texttt{srf=1.0 1.2 1.3 1.4}.
\subsection{Line types}
@@ -529,13 +520,13 @@
}}
\\
-There are only two additional informations, i.e., number of processors \texttt{\sq{nproc}} in line \texttt{\sq{preinfo}} and file name for ghost partition interfaces \texttt{\sq{gfile}} in line \texttt{\sq{mesh}} in parallel input file.
+There are only two additional pieces of information, i.e., number of processors \texttt{\sq{nproc}} in line \texttt{\sq{preinfo}} and file name for ghost partition interfaces \texttt{\sq{gfile}} in line \texttt{\sq{mesh}} in parallel input file.
\section{Input files detail}
All local element/face/edge/node numbering follows the EXODUS II convention.\\
\subsection{Coordinates files: \texttt{xfile, yfile, zfile}}
-Each of the coordinates files contains list of corresponding coordinates in following format:\\
+Each of the coordinates files contains a list of corresponding coordinates in the following format:\\
\emph{number of points \\
coordinate of point 1\\
@@ -560,7 +551,7 @@
\subsection{Connectivity file: \texttt{confile}}
-The connectivity file contains the connectivity lists of elements in following format:\\
+The connectivity file contains the connectivity lists of elements in the following format:\\
\emph{number of elements\\
$n_1$ $n_2$ $n_3$ $n_4$ $n_5$ $n_6$ $n_7$ $n_8$ of element 1\\
@@ -592,7 +583,7 @@
\subsection{Element IDs (or Material IDs) file: \texttt{idfile}}
-This file contains the IDs of elements. This ID will be used in the program mainly to identify the material regions. This file has a following format: \\
+This file contains the IDs of elements. This ID will be used in the program mainly to identify the material regions. This file has the following format: \\
\emph{number of elements\\
ID of element 1\\
@@ -625,7 +616,7 @@
\subsection{Displacement boundary conditions files: \texttt{uxfile, uyfile, uzfile}}
-This file contains information on the displacement boundary conditions (currently only the zero-displacement is implemented), and has following format:\\
+This file contains information on the displacement boundary conditions (currently only the zero-displacement is implemented), and has the following format:\\
\emph{number of element faces}\\
\emph{elementID faceID \\
@@ -665,12 +656,12 @@
...\\
...\\}
-This can be repeated as many times as many tractions.\\
+This can be repeated as many times as there are tractions.\\
-The \emph{relevant-axis} denotes the axis along which the load is varying, and it is represented by an integer as 1 = $x$-axis, 2 = $y$-axis, and 3 = $z$-axis. The variables $x_1$ and $x_2$ denote the coordinates (only the \emph{relevant-axis}) of two points between which the linearly distributed load is applied. Similarly, $q_{x1}$, $q_{y1}$ and $q_{z1}$, and $q_{x2}$, $q_{y2}$ and $q_{z2}$ denote the load vectors in kN/m\tsup{2} at the point 1 and 2, respectively.\\
+The \emph{relevant-axis} denotes the axis along which the load is varying, and is represented by an integer as 1 = $x$-axis, 2 = $y$-axis, and 3 = $z$-axis. The variables $x_1$ and $x_2$ denote the coordinates (only the \emph{relevant-axis}) of two points between which the linearly distributed load is applied. Similarly, $q_{x1}$, $q_{y1}$ and $q_{z1}$, and $q_{x2}$, $q_{y2}$ and $q_{z2}$ denote the load vectors in kN/m\tsup{2} at the point 1 and 2, respectively.\\
Example:\\
-Following data specify the two tractions: a uniformly distributed traction and a linearly distributed traction.\\\\
+The following data specify the two tractions: a uniformly distributed traction and a linearly distributed traction.\\\\
\texttt{1\\
0.0 0.0 -167.751\\
@@ -701,7 +692,7 @@
\subsection{Material list file: \texttt{matfile}}
-This file contains material properties of each material regions. Material properties must be listed in a sequential order of the unique material IDs. In addition, this data file optionally contains the information on the water condition of material regions. Material regions or material IDs must be consistent with the Material IDs (Element IDs) defined in \texttt{idfile}. The \texttt{matfile} has following format:\\\\
+This file contains material properties of each material regions. Material properties must be listed in a sequential order of the unique material IDs. In addition, this data file optionally contains the information on the water condition of material regions. Material regions or material IDs must be consistent with the Material IDs (Element IDs) defined in \texttt{idfile}. The \texttt{matfile} has the following format:\\\\
\emph{comment line}\\
\emph{number of material regions (unique material IDs)\\
@@ -720,7 +711,7 @@
Example:\\
\\
-Following data defines four material regions. No region is submerged in water.\\
+The following data defines four material regions. No region is submerged in water.\\
\texttt{\# material properties (id, domain, gamma, ym, nu, phi, coh, psi)\\
4\\
@@ -730,7 +721,7 @@
4 1 18.5 1e5 0.3 20.0 29.0 0.0\\
}\\
-Following data defines four material regions with two of them submerged.\\
+The following data defines four material regions with two of them submerged.\\
\texttt{\# material properties (id, domain, gamma, ym, nu, phi, coh, psi)\\
4\\
@@ -779,15 +770,15 @@
\subsection{Summary file}
-This file is self explanatory and it contains the summary of the result including control parameters, maximum displacement at each step, and elapsed time. The file is written in ASCII format and its name follows the convention \emph{input\_file\_name\_header}\texttt{\_summary} for serial run and \emph{input\_file\_name\_header}\texttt{\_summary\_proc}\emph{processor\_ID} for parallel run.
+This file is self explanatory and it contains a summary of the results including control parameters, maximum displacement at each step, and elapsed time. The file is written in ASCII format and its name follows the convention \emph{input\_file\_name\_header}\texttt{\_summary} for serial run and \emph{input\_file\_name\_header}\texttt{\_summary\_proc}\emph{processor\_ID} for parallel run.
\subsection{Mesh files}
-This file contains the mesh information of the model including coordinates, connectivity, element types etc. in EnSight Gold binary format~\citep[see][]{ensight2008}. The file name follows the format \emph{input\_file\_name\_header}\texttt{\_summary} for serial run and \emph{input\_file\_name\_header}\texttt{\_summary\_proc}\emph{processor\_ID} for parallel run.
+This file contains the mesh information of the model including coordinates, connectivity, element types, etc., in EnSight Gold binary format~\citep[see][]{ensight2008}. The file name follows the format \emph{input\_file\_name\_header}\texttt{\_summary} for serial run and \emph{input\_file\_name\_header}\texttt{\_summary\_proc}\emph{processor\_ID} for parallel run.
\subsection{Displacement field file}
-This file contains the nodal displacement field in the model written in EnSight Gold binary format. The file name follows the format \emph{input\_file\_name\_header}\texttt{\_step}\emph{step}\texttt{.dis} for serial run and \emph{input\_file\_name\_header}\texttt{\_step}\emph{step}\texttt{\_proc}\emph{processor\_ID}\texttt{.dis} for parallel run.
+This file contains the nodal displacement field in the model written in EnSight Gold binary format. The file name follows the format \emph{input\_file\_name\_header}\texttt{\_step}\emph{step}\texttt{.dis} for serial run and \emph{input\_file\_name\_header}\texttt{\_step}\emph{step}\texttt{\_proc}\emph{processor\_ID}\texttt{.dis} for parallel runs.
\subsection{Pore pressure file}
@@ -846,18 +837,18 @@
\texttt{gcc -o exodus2sem exodus2sem.c}
\subsubsection*{Run}
\texttt{exodus2sem} {\emph{EXODUS\_mesh\_file}} {\emph{OPTIONS}}\\
- For more detail, see \texttt{/utilities/README\_exodus2sem}. It can also be compiled automatically during the build process of main package \pack\ (see Section~\ref{sec:configure}).
+ For more details, see \texttt{/utilities/README\_exodus2sem}. It can also be compiled automatically during the build process of main package \pack\ (see Section~\ref{sec:configure}).
\section{Generate SOS file}
\label{sec:sos}
-The program \texttt{write\_sos.f90} contained in the utilities directory can be used to write EnSight Gold server-of-server file (.sos file, see~\citep{ensight2008}) to visualize the multi-processors data in parallel. This file does not contain the actual data, but only the information on the data location and parallel processing.
+The program \texttt{write\_sos.f90} contained in the utilities directory can be used to write EnSight Gold server-of-server file (.sos file, see~\citep{ensight2008}) to visualize the multi-processors data in parallel. This file does not contain the actual data, but only information on the data location and parallel processing.
\subsubsection*{Compile}
\texttt{gfortran -o write\_sos write\_sos.f90}
\subsubsection*{Run}
\texttt{exodus2sem} {\it{input\_file}}\\
-For more detail, see \texttt{/utilities/README\_write\_sos}. It can also be compiled automatically during the build process of main package \pack\ (see Section~\ref{sec:configure}).
+For more details, see \texttt{/utilities/README\_write\_sos}. It can also be compiled automatically during the build process of main package \pack\ (see Section~\ref{sec:configure}).
\bibliographystyle{chicago}
\bibliography{manual_SPECFEM3D_GEOTECH}
Modified: seismo/3D/SPECFEM3D_GEOTECH/trunk/src/Makefile
===================================================================
--- seismo/3D/SPECFEM3D_GEOTECH/trunk/src/Makefile 2012-04-21 00:26:41 UTC (rev 19963)
+++ seismo/3D/SPECFEM3D_GEOTECH/trunk/src/Makefile 2012-04-23 16:23:33 UTC (rev 19964)
@@ -9,7 +9,7 @@
#FFLAGS = -g -std=f2003 -fbounds-check -fimplicit-none -O5 -pedantic -pedantic-errors -Wline-truncation -fno-trapping-math -fno-backslash
# parallel fortran compiler
-FC = ~/openmpi-1.4.3/bin/mpif90 # OR path to appropriate MPI forrtan compiler
+FC = mpif90 # OR path to appropriate MPI forrtan compiler
FFLAGS = -g -std=f2003 -fbounds-check -fimplicit-none -O5 -pedantic -pedantic-errors -Wline-truncation -fno-trapping-math -fno-backslash
FC_COMPILE = $(FC) $(FFLAGS)
Modified: seismo/3D/SPECFEM3D_GEOTECH/trunk/src/apply_traction.f90
===================================================================
--- seismo/3D/SPECFEM3D_GEOTECH/trunk/src/apply_traction.f90 2012-04-21 00:26:41 UTC (rev 19963)
+++ seismo/3D/SPECFEM3D_GEOTECH/trunk/src/apply_traction.f90 2012-04-23 16:23:33 UTC (rev 19964)
@@ -1,3 +1,4 @@
+! TODO: make use of the orthogonality
! this routine read and applies the traction specified in the traction file
! REVISION
! HNG, Jul 12,2011; HNG, Apr 09,2010; HNG, Dec 08,2010
Modified: seismo/3D/SPECFEM3D_GEOTECH/trunk/src/compute_bmat.f90
===================================================================
--- seismo/3D/SPECFEM3D_GEOTECH/trunk/src/compute_bmat.f90 2012-04-21 00:26:41 UTC (rev 19963)
+++ seismo/3D/SPECFEM3D_GEOTECH/trunk/src/compute_bmat.f90 2012-04-23 16:23:33 UTC (rev 19964)
@@ -1,5 +1,4 @@
! This subroutine forms the strain-displacement matrix (bmat)
-! in 2D (ih=3 or 4) or 3D (ih=6)
! REFERENCE:
! copied and modified from
! Smith and Griffiths (2004): Programming the finite element method
@@ -7,47 +6,38 @@
! HNG, Jul 12,2011; HNG, Apr 09,2010
SUBROUTINE compute_bmat(bmat,deriv)
use set_precision
- IMPLICIT NONE
- REAL(kind=kreal),INTENT(IN)::deriv(:,:)
- REAL(kind=kreal),INTENT(OUT)::bmat(:,:)
- INTEGER::k,l,m,n,ih,nod
- REAL(kind=kreal)::x,y,z
- bmat=0.0_kreal
- ih=UBOUND(bmat,1)
- nod=UBOUND(deriv,2)
- SELECT CASE (ih)
- CASE(3,4)
- DO m=1,nod
- k=2*m
- l=k-1
- x=deriv(1,m)
- y=deriv(2,m)
- bmat(1,l)=x
- bmat(3,k)=x
- bmat(2,k)=y
- bmat(3,l)=y
- END DO
- CASE(6)
- DO m=1,nod
- n=3*m
- k=n-1
- l=k-1
- x=deriv(1,m)
- y=deriv(2,m)
- z=deriv(3,m)
- bmat(1,l)=x
- bmat(4,k)=x
- bmat(6,n)=x
- bmat(2,k)=y
- bmat(4,l)=y
- bmat(5,n)=y
- bmat(3,n)=z
- bmat(5,k)=z
- bmat(6,l)=z
- END DO
- CASE DEFAULT
- WRITE(*,*)'ERROR: wrong dimension for "nst" in bmat matrix!'
- END SELECT
+implicit none
+real(kind=kreal),intent(in) :: deriv(:,:)
+real(kind=kreal),intent(OUT) :: bmat(:,:)
+integer :: k,l,m,n,nod,nst
+real(kind=kreal) :: x,y,z
+
+! check size
+nst=ubound(bmat,1)
+if(nst.ne.6)then
+ write(*,*)'ERROR: wrong size of the stress tensor!'
+ stop
+endif
+nod=ubound(deriv,2)
+
+bmat=0.0_kreal
+DO m=1,nod
+ n=3*m
+ k=n-1
+ l=k-1
+ x=deriv(1,m)
+ y=deriv(2,m)
+ z=deriv(3,m)
+ bmat(1,l)=x
+ bmat(4,k)=x
+ bmat(6,n)=x
+ bmat(2,k)=y
+ bmat(4,l)=y
+ bmat(5,n)=y
+ bmat(3,n)=z
+ bmat(5,k)=z
+ bmat(6,l)=z
+END DO
RETURN
END SUBROUTINE compute_bmat
Modified: seismo/3D/SPECFEM3D_GEOTECH/trunk/src/compute_cmat.f90
===================================================================
--- seismo/3D/SPECFEM3D_GEOTECH/trunk/src/compute_cmat.f90 2012-04-21 00:26:41 UTC (rev 19963)
+++ seismo/3D/SPECFEM3D_GEOTECH/trunk/src/compute_cmat.f90 2012-04-23 16:23:33 UTC (rev 19964)
@@ -1,5 +1,4 @@
-! this subroutine returns the elastic matrix for ih=3 (plane strain),
-! ih=4 (axisymmetry or plane strain elastoplasticity) or ih=6 (3D)
+! this subroutine returns the elastic matrix (3D)
! REFERENCE:
! copied and modified from
! Smith and Griffiths (2004): Programming the finite element method
@@ -7,51 +6,37 @@
! HNG, Jul 12,2011; HNG, Apr 09,2010
subroutine compute_cmat(cmat,e,v)
use set_precision
- implicit none
- real(kind=kreal),intent(in)::e,v
- real(kind=kreal),intent(out)::cmat(:,:)
- real(kind=kreal)::v1,v2,c,vv,zero=0.0_kreal,pt5=0.5_kreal,one=1.0_kreal,two=2.0_kreal
- integer::i,ih
- cmat=zero
- ih=ubound(cmat,1)
- v1=one-v
- c=e/((one+v)*(one-two*v))
- select case(ih)
- case(3)
- cmat(1,1)=v1*c
- cmat(2,2)=v1*c
- cmat(1,2)=v*c
- cmat(2,1)=v*c
- cmat(3,3)=pt5*c*(one-two*v)
- case(4)
- cmat(1,1)=v1*c
- cmat(2,2)=v1*c
- cmat(4,4)=v1*c
- cmat(3,3)=pt5*c*(one-two*v)
- cmat(1,2)=v*c
- cmat(2,1)=v*c
- cmat(1,4)=v*c
- cmat(4,1)=v*c
- cmat(2,4)=v*c
- cmat(4,2)=v*c
- case(6)
- v2=v/(one-v)
- vv=(one-two*v)/(one-v)*pt5
- do i=1,3
- cmat(i,i)=one
- end do
- do i=4,6
- cmat(i,i)=vv
- end do
- cmat(1,2)=v2
- cmat(2,1)=v2
- cmat(1,3)=v2
- cmat(3,1)=v2
- cmat(2,3)=v2
- cmat(3,2)=v2
- cmat=cmat*e/(two*(one+v)*vv)
- case default
- write(*,*)'ERROR: wrong size for "cmat" matrix!'
- end select
+implicit none
+real(kind=kreal),intent(in)::e,v
+real(kind=kreal),intent(out)::cmat(:,:)
+real(kind=kreal)::v1,v2,c,vv,zero=0.0_kreal,pt5=0.5_kreal,one=1.0_kreal,two=2.0_kreal
+integer::i,nst
+
+! check size
+nst=ubound(cmat,1)
+if(nst.ne.6)then
+ write(*,*)'ERROR: wrong size of the stress tensor!'
+ stop
+endif
+
+cmat=zero
+v1=one-v
+c=e/((one+v)*(one-two*v))
+
+v2=v/(one-v)
+vv=(one-two*v)/(one-v)*pt5
+do i=1,3
+ cmat(i,i)=one
+end do
+do i=4,6
+ cmat(i,i)=vv
+end do
+cmat(1,2)=v2
+cmat(2,1)=v2
+cmat(1,3)=v2
+cmat(3,1)=v2
+cmat(2,3)=v2
+cmat(3,2)=v2
+cmat=cmat*e/(two*(one+v)*vv)
return
end subroutine compute_cmat
Modified: seismo/3D/SPECFEM3D_GEOTECH/trunk/src/excavation.f90
===================================================================
--- seismo/3D/SPECFEM3D_GEOTECH/trunk/src/excavation.f90 2012-04-21 00:26:41 UTC (rev 19963)
+++ seismo/3D/SPECFEM3D_GEOTECH/trunk/src/excavation.f90 2012-04-23 16:23:33 UTC (rev 19964)
@@ -3,7 +3,7 @@
! REVISION:
! April 09,2010, Hom Nath Gharti
! FEEDBACK:
-! homnath_AT_norsar_DOT_no
+! hgharti_AT_princeton_DOT_edu
module excavation
use set_precision
contains
@@ -24,7 +24,6 @@
ismat_off=.false.
ismat_off(excavid)=.true.
-
! find intact and void elements
ielmt_intact=0; ielmt_void=0
do i_elmt=1,nelmt
@@ -145,7 +144,7 @@
real(kreal) :: coord(ngnod,ndim),jac(ndim,ndim),deriv(ndim,nenod), &
bmat(nst,nedof),eld(nedof),bload(nedof),eload(nedof),sigma(nst)
integer :: egdof(nedof),num(nenod)
-integer :: i,i_elmt
+integer :: i,idof,i_elmt
! compute excavation load
do i_elmt=1,nelmt
@@ -153,7 +152,7 @@
num=g_num(:,i_elmt)
coord=transpose(g_coord(:,num(gnod))) !transpose(g_coord(:,num(1:ngnod)))
egdof=gdof_elmt(:,i_elmt) !reshape(gdof(:,g_num(:,ielmt)),(/nndof*nenod/))
-
+ idof=0
do i=1,ngll
!call shape_function(fun,gll_points(i))
! compute Jacobian at GLL point using 20 noded element
@@ -166,7 +165,10 @@
sigma=stress_local(:,i,i_elmt)
eload=MATMUL(sigma,bmat)
bload=bload+eload*detjac*gll_weights(i)
- eld(3:nedof:3)=eld(3:nedof:3)+lagrange_gll(i,:)*detjac*gll_weights(i)
+ ! since interpolation funtions are orthogonal we compute only nonzero terms
+ idof=idof+3
+ eld(idof)=eld(idof)+detjac*gll_weights(i)
+ !eld(3:nedof:3)=eld(3:nedof:3)+lagrange_gll(i,:)*detjac*gll_weights(i)
end do ! i=1,ngll
extload(egdof)=extload(egdof)+eld*gam(mat_id(i_elmt))+bload
enddo
@@ -194,7 +196,7 @@
real(kreal) :: coord(ngnod,ndim),jac(ndim,ndim),deriv(ndim,nenod), &
bmat(nst,nedof),eld(nedof),bload(nedof),eload(nedof),tload(nedof),sigma(nst)
integer :: egdof(nedof),num(nenod)
-integer :: i,i_elmt
+integer :: i,idof,i_elmt
excavload=zero
! compute excavation load
@@ -203,7 +205,7 @@
num=g_num(:,i_elmt)
coord=transpose(g_coord(:,num(gnod))) !transpose(g_coord(:,num(1:ngnod)))
!egdof=gdof_elmt(:,i_elmt) !reshape(gdof(:,g_num(:,ielmt)),(/nndof*nenod/))
-
+ idof=0
do i=1,ngll
!call shape_function(fun,gll_points(i))
! compute Jacobian at GLL point using 20 noded element
@@ -216,7 +218,10 @@
sigma=stress_local(:,i,i_elmt)
eload=MATMUL(sigma,bmat)
bload=bload+eload*detjac*gll_weights(i)
- eld(3:nedof:3)=eld(3:nedof:3)+lagrange_gll(i,:)*detjac*gll_weights(i)
+ ! since interpolation funtions are orthogonal we compute only nonzero terms
+ idof=idof+3
+ eld(idof)=eld(idof)+detjac*gll_weights(i)
+ !eld(3:nedof:3)=eld(3:nedof:3)+lagrange_gll(i,:)*detjac*gll_weights(i)
end do ! i=1,ngll
tload=eld*gam(mat_id(i_elmt))+bload
!excavload(:,num)=excavload(:,num)+reshape(tload,(/nndof,ngll/))
Modified: seismo/3D/SPECFEM3D_GEOTECH/trunk/src/gll_library.f90
===================================================================
--- seismo/3D/SPECFEM3D_GEOTECH/trunk/src/gll_library.f90 2012-04-21 00:26:41 UTC (rev 19963)
+++ seismo/3D/SPECFEM3D_GEOTECH/trunk/src/gll_library.f90 2012-04-23 16:23:33 UTC (rev 19964)
@@ -56,11 +56,11 @@
eta=gll_points(2,ii)
zeta=gll_points(3,ii)
- ! compute 1d lagrange polynomials
- call lagrange1d(ngllx,xi,lagrange_x,lagrange_dx)
- call lagrange1d(nglly,eta,lagrange_y,lagrange_dy)
- call lagrange1d(ngllz,zeta,lagrange_z,lagrange_dz)
-
+ ! compute 1d lagrange polynomials on GLL points
+ ! this can also be computed in a simple manner due to the orthogonality
+ call lagrange1dGLL(ngllx,gllpx,xi,lagrange_x,lagrange_dx)
+ call lagrange1dGLL(nglly,gllpy,eta,lagrange_y,lagrange_dy)
+ call lagrange1dGLL(ngllz,gllpz,zeta,lagrange_z,lagrange_dz)
n=0
do k=1,ngllz
do j=1,nglly
@@ -74,7 +74,6 @@
enddo
enddo
enddo
-
return
end subroutine gll_quadrature
!===========================================
@@ -127,9 +126,9 @@
eta=gll_points2d(2,ii)
! compute 1d lagrange polynomials
- call lagrange1d(ngllx,xi,lagrange_x,lagrange_dx)
- call lagrange1d(nglly,eta,lagrange_y,lagrange_dy)
-
+ ! this can also be computed in a simple manner due to the orthogonality
+ call lagrange1dGLL(ngllx,gllpx,xi,lagrange_x,lagrange_dx)
+ call lagrange1dGLL(nglly,gllpy,eta,lagrange_y,lagrange_dy)
n=0
do j=1,nglly
do i=1,ngllx
@@ -186,7 +185,8 @@
xi=gll_points1d(1,ii)
! compute 1d lagrange polynomials
- call lagrange1d(ngllx,xi,lagrange_x,lagrange_dx)
+ ! this can also be computed in a simple manner due to the orthogonality
+ call lagrange1dGLL(ngllx,gllpx,xi,lagrange_x,lagrange_dx)
n=0
do i=1,ngllx
@@ -202,7 +202,7 @@
! this subroutine computes the 1d lagrange interpolation functions and their
! derivatives at a given point xi.
-subroutine lagrange1d(nenod,xi,phi,dphi_dxi)
+subroutine lagrange1dGEN(nenod,xi,phi,dphi_dxi)
implicit none
integer,intent(in) :: nenod ! number of nodes in an 1d element
integer :: i,j,k
@@ -252,10 +252,57 @@
enddo
return
-end subroutine lagrange1d
+end subroutine lagrange1dGEN
!===========================================
+! this subroutine computes the 1d lagrange interpolation functions and their
+! derivatives at a given point xi.
+subroutine lagrange1dGLL(nenod,xii,xi,phi,dphi_dxi)
+implicit none
+integer,intent(in) :: nenod ! number of nodes in an 1d element
+real(kind=kreal),dimension(nenod),intent(in) :: xii
+real(kind=kreal),intent(in) :: xi ! point where to calculate lagrange function and
+!its derivative
+real(kind=kreal),dimension(nenod),intent(out) :: phi,dphi_dxi
+
+integer :: i,j,k
+real(kind=kreal),dimension(nenod) :: term,dterm,sum_term
+real(kind=kreal) :: dx
+
+do i=1,nenod
+ k=0
+ phi(i)=1.0_kreal
+ do j=1,nenod
+ if(j/=i)then
+ k=k+1
+ term(k)=(xi-xii(j))/(xii(i)-xii(j))
+ dterm(k)=1.0_kreal/(xii(i)-xii(j)) ! derivative of the term wrt xi
+
+ phi(i)=phi(i)*(xi-xii(j))/(xii(i)-xii(j))
+ endif
+ enddo
+
+ sum_term=1.0_kreal
+ do j=1,nenod-1
+ do k=1,nenod-1
+ if(k==j)then
+ sum_term(j)=sum_term(j)*dterm(k)
+ else
+ sum_term(j)=sum_term(j)*term(k)
+ endif
+ enddo
+ enddo
+ dphi_dxi(i)=0.0_kreal
+ do j=1,nenod-1
+ dphi_dxi(i)=dphi_dxi(i)+sum_term(j)
+ enddo
+enddo
+
+return
+end subroutine lagrange1dGLL
!===========================================
+
+!===========================================
!
! Library to compute the Gauss-Lobatto-Legendre points and weights
! Based on Gauss-Lobatto routines from M.I.T.
Modified: seismo/3D/SPECFEM3D_GEOTECH/trunk/src/global.f90
===================================================================
--- seismo/3D/SPECFEM3D_GEOTECH/trunk/src/global.f90 2012-04-21 00:26:41 UTC (rev 19963)
+++ seismo/3D/SPECFEM3D_GEOTECH/trunk/src/global.f90 2012-04-23 16:23:33 UTC (rev 19964)
@@ -46,9 +46,9 @@
integer,parameter :: nst=6 ! number of unique stress components
! sx,sy,sz,tauxy,tauyz,tauzx
-character(len=150) :: file_head,inp_path,out_path,part_path
+character(len=250) :: file_head,inp_path,out_path,part_path
! displacement BC, ghost, traction, and water surface files
-character(len=150) :: uxfile,uyfile,uzfile,gfile,trfile,wsfile
+character(len=250) :: uxfile,uyfile,uzfile,gfile,trfile,wsfile
integer :: cg_maxiter,nl_maxiter,nexcav,ninc,nsrf,ntstep
real(kind=kreal) :: cg_tol,nl_tol
integer,allocatable :: excavid(:),nexcavid(:) ! Excavation ID (regions), nunber of excavation IDs (regions) in each stage
Modified: seismo/3D/SPECFEM3D_GEOTECH/trunk/src/math_library.f90
===================================================================
--- seismo/3D/SPECFEM3D_GEOTECH/trunk/src/math_library.f90 2012-04-21 00:26:41 UTC (rev 19963)
+++ seismo/3D/SPECFEM3D_GEOTECH/trunk/src/math_library.f90 2012-04-23 16:23:33 UTC (rev 19964)
@@ -76,26 +76,26 @@
!=======================================================
! this function returns the determinant of a 1x1, 2x2 or 3x3
-! jacobian matrix.
+! matrix.
! this routine was copied and modified from
! Smith and Griffiths (2004): Programming the finite element method
-function determinant(jac)result(det)
+function determinant(xmat)result(det)
implicit none
-real(kind=kreal),intent(in)::jac(:,:)
+real(kind=kreal),intent(in) :: xmat(:,:)
real(kind=kreal)::det
-integer::it
-it=ubound(jac,1)
-select case(it)
+integer::n
+n=ubound(xmat,1)
+select case(n)
case(1)
det=1.0_kreal
case(2)
- det=jac(1,1)*jac(2,2)-jac(1,2)*jac(2,1)
+ det=xmat(1,1)*xmat(2,2)-xmat(1,2)*xmat(2,1)
case(3)
- det=jac(1,1)*(jac(2,2)*jac(3,3)-jac(3,2)*jac(2,3))
- det=det-jac(1,2)*(jac(2,1)*jac(3,3)-jac(3,1)*jac(2,3))
- det=det+jac(1,3)*(jac(2,1)*jac(3,2)-jac(3,1)*jac(2,2))
+ det=xmat(1,1)*(xmat(2,2)*xmat(3,3)-xmat(3,2)*xmat(2,3))
+ det=det-xmat(1,2)*(xmat(2,1)*xmat(3,3)-xmat(3,1)*xmat(2,3))
+ det=det+xmat(1,3)*(xmat(2,1)*xmat(3,2)-xmat(3,1)*xmat(2,2))
case default
- write(*,*)'ERROR: wrong dimension for jacobian matrix!'
+ write(*,*)'ERROR: determinant of ',n,'X',n,' matrix not supported!'
end select
return
end function determinant
@@ -104,53 +104,48 @@
! this subroutine inverts a small square matrix onto itself.
! this routine was copied and modified from
! Smith and Griffiths (2004): Programming the finite element method
-subroutine invert(matrix)
+subroutine invert(xmat)
implicit none
-real(kind=kreal),intent(in out)::matrix(:,:)
+real(kind=kreal),intent(in out)::xmat(:,:)
real(kind=kreal)::det,j11,j12,j13,j21,j22,j23,j31,j32,j33,con
integer::ndim,i,k
-ndim=ubound(matrix,1)
+ndim=ubound(xmat,1)
if(ndim==2)then
- det=matrix(1,1)*matrix(2,2)-matrix(1,2)*matrix(2,1)
- j11=matrix(1,1)
- matrix(1,1)=matrix(2,2)
- matrix(2,2)=j11
- matrix(1,2)=-matrix(1,2)
- matrix(2,1)=-matrix(2,1)
- matrix=matrix/det
+ det=xmat(1,1)*xmat(2,2)-xmat(1,2)*xmat(2,1)
+ j11=xmat(1,1)
+ xmat(1,1)=xmat(2,2)
+ xmat(2,2)=j11
+ xmat(1,2)=-xmat(1,2)
+ xmat(2,1)=-xmat(2,1)
+ xmat=xmat/det
else if(ndim==3)then
- det=matrix(1,1)*(matrix(2,2)*matrix(3,3)-matrix(3,2)*matrix(2,3))
- det=det-matrix(1,2)*(matrix(2,1)*matrix(3,3)-matrix(3,1)*matrix(2,3))
- det=det+matrix(1,3)*(matrix(2,1)*matrix(3,2)-matrix(3,1)*matrix(2,2))
- j11=matrix(2,2)*matrix(3,3)-matrix(3,2)*matrix(2,3)
- j21=-matrix(2,1)*matrix(3,3)+matrix(3,1)*matrix(2,3)
- j31=matrix(2,1)*matrix(3,2)-matrix(3,1)*matrix(2,2)
- j12=-matrix(1,2)*matrix(3,3)+matrix(3,2)*matrix(1,3)
- j22=matrix(1,1)*matrix(3,3)-matrix(3,1)*matrix(1,3)
- j32=-matrix(1,1)*matrix(3,2)+matrix(3,1)*matrix(1,2)
- j13=matrix(1,2)*matrix(2,3)-matrix(2,2)*matrix(1,3)
- j23=-matrix(1,1)*matrix(2,3)+matrix(2,1)*matrix(1,3)
- j33=matrix(1,1)*matrix(2,2)-matrix(2,1)*matrix(1,2)
- matrix(1,1)=j11
- matrix(1,2)=j12
- matrix(1,3)=j13
- matrix(2,1)=j21
- matrix(2,2)=j22
- matrix(2,3)=j23
- matrix(3,1)=j31
- matrix(3,2)=j32
- matrix(3,3)=j33
- matrix=matrix/det
+ det=xmat(1,1)*(xmat(2,2)*xmat(3,3)-xmat(3,2)*xmat(2,3))
+ det=det-xmat(1,2)*(xmat(2,1)*xmat(3,3)-xmat(3,1)*xmat(2,3))
+ det=det+xmat(1,3)*(xmat(2,1)*xmat(3,2)-xmat(3,1)*xmat(2,2))
+ j11=xmat(2,2)*xmat(3,3)-xmat(3,2)*xmat(2,3)
+ j21=-xmat(2,1)*xmat(3,3)+xmat(3,1)*xmat(2,3)
+ j31=xmat(2,1)*xmat(3,2)-xmat(3,1)*xmat(2,2)
+ j12=-xmat(1,2)*xmat(3,3)+xmat(3,2)*xmat(1,3)
+ j22=xmat(1,1)*xmat(3,3)-xmat(3,1)*xmat(1,3)
+ j32=-xmat(1,1)*xmat(3,2)+xmat(3,1)*xmat(1,2)
+ j13=xmat(1,2)*xmat(2,3)-xmat(2,2)*xmat(1,3)
+ j23=-xmat(1,1)*xmat(2,3)+xmat(2,1)*xmat(1,3)
+ j33=xmat(1,1)*xmat(2,2)-xmat(2,1)*xmat(1,2)
+
+ xmat(1,1)=j11; xmat(1,2)=j12; xmat(1,3)=j13
+ xmat(2,1)=j21; xmat(2,2)=j22; xmat(2,3)=j23
+ xmat(3,1)=j31; xmat(3,2)=j32; xmat(3,3)=j33
+ xmat=xmat/det
else
do k=1,ndim
- con=matrix(k,k)
- matrix(k,k)=1.0_kreal
- matrix(k,:)=matrix(k,:)/con
+ con=xmat(k,k)
+ xmat(k,k)=1.0_kreal
+ xmat(k,:)=xmat(k,:)/con
do i=1,ndim
if(i/=k)then
- con=matrix(i,k)
- matrix(i,k)=0.0_kreal
- matrix(i,:)=matrix(i,:)-matrix(k,:)*con
+ con=xmat(i,k)
+ xmat(i,k)=0.0_kreal
+ xmat(i,:)=xmat(i,:)-xmat(k,:)*con
end if
end do
end do
@@ -159,67 +154,50 @@
end subroutine invert
!=======================================================
-! this subroutine forms the stress invariants in 2- or 3-d.
+! this subroutine forms the stress invariants 3d.
! this routine was copied and modified from
! Smith and Griffiths (2004): Programming the finite element method
subroutine stress_invariant(stress,sigm,dsbar,theta)
implicit none
real(kind=kreal),intent(in)::stress(:)
real(kind=kreal),intent(out),optional::sigm,dsbar,theta
-real(kind=kreal)::sx,sy,sz,txy,dx,dy,dz,xj3,sine,s1,s2,s3,s4,s5,s6,ds1,ds2,ds3,&
- d2,d3,sq3,zero=0.0_kreal,small=1.e-12_kreal,one=1.0_kreal,two=2.0_kreal, &
- three=3.0_kreal,six=6.0_kreal,thpt5=13.5_kreal
-integer::nst
+real(kind=kreal)::sx,sy,sz,txy,dx,dy,dz,xj3,sine,s1,s2,s3,s4,s5,s6,ds1,ds2, &
+ds3,d2,d3,sq3,zero=0.0_kreal,small=1.e-12_kreal,one=1.0_kreal,two=2.0_kreal, &
+three=3.0_kreal,six=6.0_kreal,thpt5=13.5_kreal
+integer :: nst
+
+! check size
nst=ubound(stress,1)
-select case(nst)
-case(4)
- sx=stress(1)
- sy=stress(2)
- txy=stress(3)
- sz=stress(4)
- sigm=(sx+sy+sz)/three
- dsbar=sqrt((sx-sy)**2+(sy-sz)**2+(sz-sx)**2+six*txy**2)/sqrt(two)
- if(dsbar<small)then
- theta=zero
- else
- dx=(two*sx-sy-sz)/three
- dy=(two*sy-sz-sx)/three
- dz=(two*sz-sx-sy)/three
- xj3=dx*dy*dz-dz*txy**2
- sine=-thpt5*xj3/dsbar**3
- if(sine>=one)sine=one
- if(sine<-one)sine=-one
- theta=asin(sine)/three
- end if
-case(6)
- sq3=sqrt(three)
- s1=stress(1)
- s2=stress(2)
- s3=stress(3)
- s4=stress(4)
- s5=stress(5)
- s6=stress(6)
- sigm=(s1+s2+s3)/three
- d2=((s1-s2)**2+(s2-s3)**2+(s3-s1)**2)/six+s4*s4+s5*s5+s6*s6
+if(nst.ne.6)then
+ write(*,*)'ERROR: wrong size of the stress tensor!'
+ stop
+endif
+
+sq3=sqrt(three)
+s1=stress(1)
+s2=stress(2)
+s3=stress(3)
+s4=stress(4)
+s5=stress(5)
+s6=stress(6)
+sigm=(s1+s2+s3)/three
+d2=((s1-s2)**2+(s2-s3)**2+(s3-s1)**2)/six+s4*s4+s5*s5+s6*s6
- if(d2<small)d2=small ! special case of hydrostatic pressure or just at the tip
+if(d2<small)d2=small ! special case of hydrostatic pressure or just at the tip
- ds1=s1-sigm
- ds2=s2-sigm
- ds3=s3-sigm
- d3=ds1*ds2*ds3-ds1*s5*s5-ds2*s6*s6-ds3*s4*s4+two*s4*s5*s6
- dsbar=sq3*sqrt(d2)
- if(dsbar<small)then
- theta=zero
- else
- sine=-three*sq3*d3/(two*sqrt(d2)**3)
- if(sine>=one)sine=one
- if(sine<-one)sine=-one
- theta=asin(sine)/three
- end if
-case default
- write(*,*)"ERROR: wrong size for nst in invar!"
-end select
+ds1=s1-sigm
+ds2=s2-sigm
+ds3=s3-sigm
+d3=ds1*ds2*ds3-ds1*s5*s5-ds2*s6*s6-ds3*s4*s4+two*s4*s5*s6
+dsbar=sq3*sqrt(d2)
+if(dsbar<small)then
+ theta=zero
+else
+ sine=-three*sq3*d3/(two*sqrt(d2)**3)
+ if(sine>=one)sine=one
+ if(sine<-one)sine=-one
+ theta=asin(sine)/three
+end if
return
end subroutine stress_invariant
!=======================================================
Modified: seismo/3D/SPECFEM3D_GEOTECH/trunk/src/plastic_library.f90
===================================================================
--- seismo/3D/SPECFEM3D_GEOTECH/trunk/src/plastic_library.f90 2012-04-21 00:26:41 UTC (rev 19963)
+++ seismo/3D/SPECFEM3D_GEOTECH/trunk/src/plastic_library.f90 2012-04-23 16:23:33 UTC (rev 19964)
@@ -115,132 +115,95 @@
! Smith and Griffiths (2004): Programming the finite element method
subroutine mohcouq(psi,dsbar,theta,dq1,dq2,dq3)
implicit none
- real(kind=kreal),intent(in)::psi,dsbar,theta
- real(kind=kreal),intent(out)::dq1,dq2,dq3
- real(kind=kreal)::psir,snth,snps,sq3,c1,csth,cs3th,tn3th,tnth,pt49=0.49_kreal,&
- pt5=0.5_kreal,r3=3.0_kreal
+real(kind=kreal),intent(in)::psi,dsbar,theta
+real(kind=kreal),intent(out)::dq1,dq2,dq3
+real(kind=kreal)::psir,snth,snps,sq3,c1,csth,cs3th,tn3th,tnth,pt49=0.49_kreal,&
+pt5=0.5_kreal,r3=3.0_kreal
- psir=psi*deg2rad
- snth=sin(theta)
- snps=sin(psir)
- sq3=sqrt(r3)
- dq1=snps
+psir=psi*deg2rad
+snth=sin(theta)
+snps=sin(psir)
+sq3=sqrt(r3)
+dq1=snps
- if(abs(snth).gt.pt49)then
- c1=one
- if(snth.lt.zero)c1=-one
- dq2=(sq3*pt5-c1*snps*pt5/sq3)*sq3*pt5/dsbar
- dq3=zero
- else
- csth=cos(theta)
- cs3th=cos(r3*theta)
- tn3th=tan(r3*theta)
- tnth=snth/csth
- dq2=sq3*csth/dsbar*((one+tnth*tn3th)+snps*(tn3th-tnth)/sq3)*pt5
- dq3=pt5*r3*(sq3*snth+snps*csth)/(cs3th*dsbar*dsbar)
- end if
+if(abs(snth).gt.pt49)then
+ c1=one
+ if(snth.lt.zero)c1=-one
+ dq2=(sq3*pt5-c1*snps*pt5/sq3)*sq3*pt5/dsbar
+ dq3=zero
+else
+ csth=cos(theta)
+ cs3th=cos(r3*theta)
+ tn3th=tan(r3*theta)
+ tnth=snth/csth
+ dq2=sq3*csth/dsbar*((one+tnth*tn3th)+snps*(tn3th-tnth)/sq3)*pt5
+ dq3=pt5*r3*(sq3*snth+snps*csth)/(cs3th*dsbar*dsbar)
+end if
return
end subroutine mohcouq
!=======================================================
! this subroutine forms the derivatives of the invariants with respect to
-! stress in 2- or 3-d.
+! stress in 3d.
! this routine was copied and modified from
! Smith and Griffiths (2004): Programming the finite element method
subroutine formm(stress,m1,m2,m3)
- implicit none
- real(kind=kreal),intent(in)::stress(:)
- real(kind=kreal),intent(out)::m1(:,:),m2(:,:),m3(:,:)
- real(kind=kreal)::sx,sy,txy,tyz,tzx,sz,dx,dy,dz,sigm, &
- r3=3.0_kreal,r6=6.0_kreal,r9=9.0_kreal
- integer::nst,i,j
- nst=ubound(stress,1)
- select case(nst)
- case(4)
- sx=stress(1); sy=stress(2); sz=stress(4)
- txy=stress(3)
+implicit none
+real(kind=kreal),intent(in)::stress(:)
+real(kind=kreal),intent(out)::m1(:,:),m2(:,:),m3(:,:)
+real(kind=kreal)::sx,sy,txy,tyz,tzx,sz,dx,dy,dz,sigm, &
+r3=3.0_kreal,r6=6.0_kreal,r9=9.0_kreal
+integer::nst,i,j
- dx=(two*sx-sy-sz)/r3; dy=(two*sy-sz-sx)/r3; dz=(two*sz-sx-sy)/r3
- sigm=(sx+sy+sz)/r3
- m1=zero; m2=zero; m3=zero
- m1(1,1:2)=one
- m1(2,1:2)=one
- m1(4,1:2)=one
- m1(1,4)=one
- m1(4,4)=one
- m1(2,4)=one
- m1=m1/r9/sigm
- m2(1,1)=two/r3
- m2(2,2)=two/r3
- m2(4,4)= two/r3
- m2(2,4)=-one/r3
- m2(4,2)=-one/r3
- m2(1,2)=-one/r3
- m2(2,1)=-one/r3
- m2(1,4)=-one/r3
- m2(4,1)=-one/r3
- m2(3,3)=two
- m3(3,3)=-dz
- m3(1:2,3)=txy/r3
- m3(3,1:2)=txy/r3
- m3(3,4)=-two*txy/r3
- m3(4,3)=-two*txy/r3
- m3(1,1)=dx/r3
- m3(2,4)=dx/r3
- m3(4,2)=dx/r3
- m3(2,2)=dy/r3
- m3(1,4)=dy/r3
- m3(4,1)=dy/r3
- m3(4,4)=dz/r3
- m3(1,2)=dz/r3
- m3(2,1)=dz/r3
- case(6)
- sx=stress(1); sy=stress(2); sz=stress(3)
- txy=stress(4); tyz=stress(5); tzx=stress(6)
- sigm=(sx+sy+sz)/r3
- dx=sx-sigm; dy=sy-sigm; dz=sz-sigm
- m1=zero; m2=zero
- m1(1:3,1:3)=one/(r3*sigm)
- do i=1,3
- m2(i,i)=two
- m2(i+3,i+3)=r6
- end do
- m2(1,2)=-one
- m2(1,3)=-one
- m2(2,3)=-one
- m3(1,1)=dx
- m3(1,2)=dz
- m3(1,3)=dy
- m3(1,4)=txy
- m3(1,5)=-two*tyz
- m3(1,6)=tzx
- m3(2,2)=dy
- m3(2,3)=dx
- m3(2,4)=txy
- m3(2,5)=tyz
- m3(2,6)=-two*tzx
- m3(3,3)=dz
- m3(3,4)=-two*txy
- m3(3,5)=tyz
- m3(3,6)=tzx
- m3(4,4)=-r3*dz
- m3(4,5)=r3*tzx
- m3(4,6)=r3*tyz
- m3(5,5)=-r3*dx
- m3(5,6)=r3*txy
- m3(6,6)=-r3*dy
- do i=1,6
- do j=i+1,6
- m1(j,i)=m1(i,j)
- m2(j,i)=m2(i,j)
- m3(j,i)=m3(i,j)
- end do
- end do
- m1=m1/r3; m2=m2/r3; m3=m3/r3
- case default
- write(*,*)"ERROR: nst size not recognized in formm!"
- stop
- end select
+nst=ubound(stress,1)
+if(nst.ne.6)then
+ write(*,*)'ERROR: wrong size of the stress tensor!'
+ stop
+endif
+
+sx=stress(1); sy=stress(2); sz=stress(3)
+txy=stress(4); tyz=stress(5); tzx=stress(6)
+sigm=(sx+sy+sz)/r3
+dx=sx-sigm; dy=sy-sigm; dz=sz-sigm
+
+m1=zero; m2=zero
+m1(1:3,1:3)=one/(r3*sigm)
+do i=1,3
+ m2(i,i)=two
+ m2(i+3,i+3)=r6
+end do
+m2(1,2)=-one
+m2(1,3)=-one
+m2(2,3)=-one
+m3(1,1)=dx
+m3(1,2)=dz
+m3(1,3)=dy
+m3(1,4)=txy
+m3(1,5)=-two*tyz
+m3(1,6)=tzx
+m3(2,2)=dy
+m3(2,3)=dx
+m3(2,4)=txy
+m3(2,5)=tyz
+m3(2,6)=-two*tzx
+m3(3,3)=dz
+m3(3,4)=-two*txy
+m3(3,5)=tyz
+m3(3,6)=tzx
+m3(4,4)=-r3*dz
+m3(4,5)=r3*tzx
+m3(4,6)=r3*tyz
+m3(5,5)=-r3*dx
+m3(5,6)=r3*txy
+m3(6,6)=-r3*dy
+do i=1,6
+ do j=i+1,6
+ m1(j,i)=m1(i,j)
+ m2(j,i)=m2(i,j)
+ m3(j,i)=m3(i,j)
+ end do
+end do
+m1=m1/r3; m2=m2/r3; m3=m3/r3
return
end subroutine formm
!=======================================================
Modified: seismo/3D/SPECFEM3D_GEOTECH/trunk/src/preprocess.f90
===================================================================
--- seismo/3D/SPECFEM3D_GEOTECH/trunk/src/preprocess.f90 2012-04-21 00:26:41 UTC (rev 19963)
+++ seismo/3D/SPECFEM3D_GEOTECH/trunk/src/preprocess.f90 2012-04-23 16:23:33 UTC (rev 19964)
@@ -44,7 +44,7 @@
real(kind=kreal) :: cmat(nst,nst),coord(ngnod,ndim),jac(ndim,ndim),deriv(ndim,nenod), &
bmat(nst,nedof),eld(nedof),eqload(nedof),km(nedof,nedof)
integer :: egdof(nedof),num(nenod)
-integer :: i,i_elmt,k
+integer :: i,idof,i_elmt,k
if(present(extload).and.(.not.present(gravity) .or. .not.present(pseudoeq)))then
write(*,'(/,a)')'ERROR: both "gravity" and "pseudoeq" must be defined for "extload"!'
@@ -61,6 +61,7 @@
!print*,ngllx,nglly,ngllz,ngll,nndof,nenod ,gdof(:,g_num(:,ielmt))
egdof=gdof_elmt(:,i_elmt) !reshape(gdof(:,g_num(:,ielmt)),(/nndof*nenod/)) !g=g_g(:,ielmt)
km=zero; eld=zero; eqload=zero
+ idof=0
do i=1,ngll
!call shape_function(fun,gll_points(i))
! compute Jacobian at GLL point using 20 noded element
@@ -72,8 +73,10 @@
deriv=matmul(jac,dlagrange_gll(:,i,:)) ! use der for gll
call compute_bmat(bmat,deriv) !!! gll bmat matrix
km=km+matmul(matmul(transpose(bmat),cmat),bmat)*detjac*gll_weights(i)
- eld(3:nedof:3)=eld(3:nedof:3)+lagrange_gll(i,:)*detjac*gll_weights(i)
- !eld(2:nedof-1:3)=eld(2:nedof-1:3)+fun(:)*detjac*weights(i)
+ idof=idof+3
+ ! interpolation functions are orthogonal, hence it is simple
+ eld(idof)=eld(idof)+detjac*gll_weights(i)
+ !eld(3:nedof:3)=eld(3:nedof:3)+lagrange_gll(i,:)*detjac*gll_weights(i)
end do ! i=1,ngll
storkm(:,:,i_elmt)=km
do k=1,nedof
Modified: seismo/3D/SPECFEM3D_GEOTECH/trunk/src/read_input.f90
===================================================================
--- seismo/3D/SPECFEM3D_GEOTECH/trunk/src/read_input.f90 2012-04-21 00:26:41 UTC (rev 19963)
+++ seismo/3D/SPECFEM3D_GEOTECH/trunk/src/read_input.f90 2012-04-23 16:23:33 UTC (rev 19964)
@@ -15,13 +15,13 @@
integer,intent(in) :: myid,nproc
integer,intent(out) :: errcode
character(len=250),intent(out) :: errtag
-character(len=256) :: line
-character(len=400) ::tag
-character(len=80) :: strval,token
+character(len=250) :: line
+character(len=800) ::tag
+character(len=250) :: strval,token
character(len=1) :: tmp_char
-character(len=80),dimension(50) :: args
-character(len=80) :: confile,idfile,matfile
-character(len=80),dimension(3) :: coordfile
+character(len=250),dimension(50) :: args
+character(len=250) :: confile,idfile,matfile
+character(len=250),dimension(3) :: coordfile
integer :: id,ind,ios,narg,slen
integer :: bc_stat,preinfo_stat,mesh_stat,material_stat,control_stat,eqload_stat, &
@@ -31,7 +31,7 @@
character(len=20) :: format_str,ptail
character(len=250) :: fname
-character(len=150) :: data_path,mat_path
+character(len=250) :: data_path,mat_path
integer :: ipart,nproc_inp ! partition ID
integer :: iselastic,ismatpart,istat,ival,issave,nexcavid_all
Modified: seismo/3D/SPECFEM3D_GEOTECH/trunk/src/semexcav3d.F90
===================================================================
--- seismo/3D/SPECFEM3D_GEOTECH/trunk/src/semexcav3d.F90 2012-04-21 00:26:41 UTC (rev 19963)
+++ seismo/3D/SPECFEM3D_GEOTECH/trunk/src/semexcav3d.F90 2012-04-23 16:23:33 UTC (rev 19964)
@@ -1,5 +1,7 @@
! include 'license.txt'
! this is a main routine for multistage excavation
+! this program was originally based on the book "Programming the finite element
+! method" Smith and Griffiths (2004)
! REVISION:
! HNG, Aug 25,2011; HNG, Jul 14,2011; HNG, Jul 11,2011; Apr 09,2010
subroutine semexcav3d(ismpi,myid,nproc,gnod,sum_file,ptail,format_str)
Modified: seismo/3D/SPECFEM3D_GEOTECH/trunk/src/semgeotech.F90
===================================================================
--- seismo/3D/SPECFEM3D_GEOTECH/trunk/src/semgeotech.F90 2012-04-21 00:26:41 UTC (rev 19963)
+++ seismo/3D/SPECFEM3D_GEOTECH/trunk/src/semgeotech.F90 2012-04-23 16:23:33 UTC (rev 19964)
@@ -23,7 +23,7 @@
integer :: gnod(8),map2exodus(8),ngllxy,node_hex8(8)
character(len=250) :: arg1,inp_fname,out_fname,prog
-character(len=150) :: path
+character(len=250) :: path
character(len=20), parameter :: wild_char='********************'
character(len=20) :: ensight_etype
character(len=80) :: buffer,destag ! this must be 80 characters long
@@ -65,7 +65,7 @@
call close_process()
elseif(trim(arg1)==('--version'))then
if(myid==1)then
- write(stdout,'(a)')'SPECFEM3D_GEOTECH 1.0 Beta'
+ write(stdout,'(a)')'SPECFEM3D_GEOTECH 1.1 Beta'
write(stdout,'(a)')'This is free software; see the source for copying '
write(stdout,'(a)')'conditions. There is NO warranty; not even for '
write(stdout,'(a)')'MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.'
Modified: seismo/3D/SPECFEM3D_GEOTECH/trunk/src/semslope3d.F90
===================================================================
--- seismo/3D/SPECFEM3D_GEOTECH/trunk/src/semslope3d.F90 2012-04-21 00:26:41 UTC (rev 19963)
+++ seismo/3D/SPECFEM3D_GEOTECH/trunk/src/semslope3d.F90 2012-04-23 16:23:33 UTC (rev 19964)
@@ -1,4 +1,6 @@
! this is a main routine for slope stabiliy analysis
+! this program was originally based on the book "Programming the finite element
+! method" Smith and Griffiths (2004)
! REVISION:
! HNG, Jul 14,2011; HNG, Jul 11,2011; Apr 09,2010
subroutine semslope3d(ismpi,myid,nproc,gnod,sum_file,ptail,format_str)
Modified: seismo/3D/SPECFEM3D_GEOTECH/trunk/src/string_library.f90
===================================================================
--- seismo/3D/SPECFEM3D_GEOTECH/trunk/src/string_library.f90 2012-04-21 00:26:41 UTC (rev 19963)
+++ seismo/3D/SPECFEM3D_GEOTECH/trunk/src/string_library.f90 2012-04-23 16:23:33 UTC (rev 19964)
@@ -193,11 +193,11 @@
! get string value from string list which contain a character '=' that separates
! variable name and variable vlue
-character(len=80) function get_string(vname,slist,nvar)
+character(len=250) function get_string(vname,slist,nvar)
character(len=*),intent(in) :: vname
character(len=*),dimension(*) :: slist
integer,intent(in) :: nvar
-character(len=80),dimension(2) :: args
+character(len=250),dimension(2) :: args
integer :: i,narg
do i=1,nvar
@@ -221,7 +221,7 @@
character(len=*),intent(out) :: strval
character(len=*),dimension(*) :: slist
integer,intent(in) :: nvar
-character(len=80),dimension(2) :: args
+character(len=250),dimension(2) :: args
integer :: i,narg
strval=''
@@ -245,11 +245,11 @@
function get_string_vect(vname,n,slist,nvar)
implicit none
integer,intent(in) :: n
-character(len=80),dimension(n) :: get_string_vect
+character(len=250),dimension(n) :: get_string_vect
character(len=*),intent(in) :: vname
character(len=*),dimension(*),intent(in) :: slist
integer,intent(in) :: nvar
-character(len=80),dimension(2) :: args
+character(len=250),dimension(2) :: args
integer :: i,narg
do i=1,nvar
@@ -272,7 +272,7 @@
character(len=*),intent(in) :: vname
character(len=*),dimension(*),intent(in) :: slist
integer,intent(in) :: nvar
-character(len=80),dimension(2) :: args
+character(len=250),dimension(2) :: args
integer :: i,narg
do i=1,nvar
@@ -296,7 +296,7 @@
character(len=*),dimension(*),intent(in) :: slist
integer,intent(in) :: nvar
integer,intent(out) :: ival,istat
-character(len=80),dimension(2) :: args
+character(len=250),dimension(2) :: args
integer :: i,narg
ival=0
istat=-1
@@ -325,7 +325,7 @@
character(len=*),intent(in) :: vname
character(len=*),dimension(*),intent(in) :: slist
integer,intent(in) :: nvar
-character(len=80),dimension(2) :: args
+character(len=250),dimension(2) :: args
integer :: i,ios,narg
do i=1,nvar
@@ -353,7 +353,7 @@
character(len=*),dimension(*),intent(in) :: slist
integer,intent(in) :: nvar
integer,intent(out) :: istat
-character(len=80),dimension(2) :: args
+character(len=250),dimension(2) :: args
integer :: i,ios,narg
ivect=0
istat=-1
@@ -379,7 +379,7 @@
character(len=*),intent(in) :: vname
character(len=*),dimension(*),intent(in) :: slist
integer,intent(in) :: nvar
-character(len=80),dimension(2) :: args
+character(len=250),dimension(2) :: args
integer :: i,narg
do i=1,nvar
@@ -404,7 +404,7 @@
integer,intent(in) :: nvar
integer,intent(out) :: istat
real(kind=kreal),intent(out) :: rval
-character(len=80),dimension(2) :: args
+character(len=250),dimension(2) :: args
integer :: i,narg
rval=0_kreal
istat=-1
@@ -428,7 +428,7 @@
character(len=*),intent(in) :: vname
character(len=*),dimension(*),intent(in) :: slist
integer,intent(in) :: nvar
-character(len=80),dimension(2) :: args
+character(len=250),dimension(2) :: args
integer :: i,narg
do i=1,nvar
@@ -454,7 +454,7 @@
character(len=*),intent(in) :: vname
character(len=*),dimension(*),intent(in) :: slist
integer,intent(in) :: nvar
-character(len=80),dimension(2) :: args
+character(len=250),dimension(2) :: args
integer :: i,ios,narg
do i=1,nvar
@@ -483,7 +483,7 @@
character(len=*),dimension(*),intent(in) :: slist
integer,intent(in) :: nvar
integer,intent(out) :: istat
-character(len=80),dimension(2) :: args
+character(len=250),dimension(2) :: args
integer :: i,ios,narg
rvect=0.0_kreal
istat=-1
@@ -506,7 +506,7 @@
!=====================================================
! get format string for intger
-character(len=80) function form4int(n)
+character(len=250) function form4int(n)
integer,intent(in) :: n
! format for integer n
More information about the CIG-COMMITS
mailing list