[cig-commits] r22850 - in seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT: . doc example geocubitlib geocubitlib/dev utility
carltape at geodynamics.org
carltape at geodynamics.org
Wed Sep 25 17:57:45 PDT 2013
Author: carltape
Date: 2013-09-25 17:57:45 -0700 (Wed, 25 Sep 2013)
New Revision: 22850
Added:
seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/LICENSE.txt
seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/README.md
seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/example/volume_from_irregularsurf.cfg
seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/dev/cmpl.py
seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/sea.py
Removed:
seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/doc/geocubit_manual.pdfsync
seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/example/data.zip
seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/dev/cubit2googleearth.py
seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/dev/material.py
seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/obsolete/
seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/save_fault_nodes_elements.py
seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/run_cubit2specfem3d.py
seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/test_cpml_2layers.cub
seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/utility/jq.sh
Modified:
seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/GEOCUBIT.py
seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/doc/geocubit_manual.pdf
seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/doc/geocubit_manual.tex
seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/boundary_definition.py
seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/cubit2specfem3d.py
seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/exportlib.py
seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/local_volume.py
seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/menu.py
seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/mesh_volume.py
seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/read_parameter_cfg.py
seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/start.py
seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/utilities.py
seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/volumes.py
seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/setup.py
Log:
modified CUBIT_GEOCUBIT/ to match clone from github on 25-Sept-2013
Modified: seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/GEOCUBIT.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/GEOCUBIT.py 2013-09-26 00:49:25 UTC (rev 22849)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/GEOCUBIT.py 2013-09-26 00:57:45 UTC (rev 22850)
@@ -55,7 +55,7 @@
from geocubitlib import mesh_volume
mesh_volume.mesh()
- #EXPORT
+ #MERGING and EXPORTING
if menu.collect:
from geocubitlib.exportlib import collect
try:
@@ -64,10 +64,10 @@
output='totalmesh_merged'
output=output.upper()
#
- print menu.cpuxmin,menu.cpuxmax,menu.cpuymin,menu.cpuymax,menu.cpux,menu.cpuy
- collect(menu.cpuxmin,menu.cpuxmax,menu.cpuymin,menu.cpuymax,menu.cpux,menu.cpuy,menu.cubfiles,menu.ckbound_method1,menu.ckbound_method2,menu.merge_tolerance,curverefining=menu.curverefining,outfilename=output,qlog=menu.qlog,export2SPECFEM3D=menu.export2SPECFEM3D,listblock=menu.listblock,listflag=menu.listflag,outdir=menu.SPECFEM3D_output_dir)
+ collect(menu.cpuxmin,menu.cpuxmax,menu.cpuymin,menu.cpuymax,menu.cpux,menu.cpuy,menu.cubfiles,menu.ckbound_method1,menu.ckbound_method2,menu.merge_tolerance,curverefining=menu.curverefining,outfilename=output,qlog=menu.qlog,export2SPECFEM3D=menu.export2SPECFEM3D,listblock=menu.listblock,listflag=menu.listflag,outdir=menu.SPECFEM3D_output_dir,add_sea=menu.add_sea,decimate=menu.decimate,cpml=menu.cpml,cpml_size=menu.cpml_size,top_absorbing=menu.top_absorbing,hex27=menu.hex27)
if menu.export2SPECFEM3D and not menu.collect:
from geocubitlib.exportlib import e2SEM
print menu.cubfiles
- e2SEM(files=menu.cubfiles,listblock=menu.listblock,listflag=menu.listflag,outdir=menu.SPECFEM3D_output_dir)
+ print 'hex27 ',menu.hex27
+ e2SEM(files=menu.cubfiles,listblock=menu.listblock,listflag=menu.listflag,outdir=menu.SPECFEM3D_output_dir,hex27=menu.hex27)
Added: seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/LICENSE.txt
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/LICENSE.txt (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/LICENSE.txt 2013-09-26 00:57:45 UTC (rev 22850)
@@ -0,0 +1,674 @@
+ GNU GENERAL PUBLIC LICENSE
+ Version 3, 29 June 2007
+
+ Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/>
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+ Preamble
+
+ The GNU General Public License is a free, copyleft license for
+software and other kinds of works.
+
+ The licenses for most software and other practical works are designed
+to take away your freedom to share and change the works. By contrast,
+the GNU General Public License is intended to guarantee your freedom to
+share and change all versions of a program--to make sure it remains free
+software for all its users. We, the Free Software Foundation, use the
+GNU General Public License for most of our software; it applies also to
+any other work released this way by its authors. You can apply it to
+your programs, too.
+
+ When we speak of free software, we are referring to freedom, not
+price. Our General Public Licenses are designed to make sure that you
+have the freedom to distribute copies of free software (and charge for
+them if you wish), that you receive source code or can get it if you
+want it, that you can change the software or use pieces of it in new
+free programs, and that you know you can do these things.
+
+ To protect your rights, we need to prevent others from denying you
+these rights or asking you to surrender the rights. Therefore, you have
+certain responsibilities if you distribute copies of the software, or if
+you modify it: responsibilities to respect the freedom of others.
+
+ For example, if you distribute copies of such a program, whether
+gratis or for a fee, you must pass on to the recipients the same
+freedoms that you received. You must make sure that they, too, receive
+or can get the source code. And you must show them these terms so they
+know their rights.
+
+ Developers that use the GNU GPL protect your rights with two steps:
+(1) assert copyright on the software, and (2) offer you this License
+giving you legal permission to copy, distribute and/or modify it.
+
+ For the developers' and authors' protection, the GPL clearly explains
+that there is no warranty for this free software. For both users' and
+authors' sake, the GPL requires that modified versions be marked as
+changed, so that their problems will not be attributed erroneously to
+authors of previous versions.
+
+ Some devices are designed to deny users access to install or run
+modified versions of the software inside them, although the manufacturer
+can do so. This is fundamentally incompatible with the aim of
+protecting users' freedom to change the software. The systematic
+pattern of such abuse occurs in the area of products for individuals to
+use, which is precisely where it is most unacceptable. Therefore, we
+have designed this version of the GPL to prohibit the practice for those
+products. If such problems arise substantially in other domains, we
+stand ready to extend this provision to those domains in future versions
+of the GPL, as needed to protect the freedom of users.
+
+ Finally, every program is threatened constantly by software patents.
+States should not allow patents to restrict development and use of
+software on general-purpose computers, but in those that do, we wish to
+avoid the special danger that patents applied to a free program could
+make it effectively proprietary. To prevent this, the GPL assures that
+patents cannot be used to render the program non-free.
+
+ The precise terms and conditions for copying, distribution and
+modification follow.
+
+ TERMS AND CONDITIONS
+
+ 0. Definitions.
+
+ "This License" refers to version 3 of the GNU General Public License.
+
+ "Copyright" also means copyright-like laws that apply to other kinds of
+works, such as semiconductor masks.
+
+ "The Program" refers to any copyrightable work licensed under this
+License. Each licensee is addressed as "you". "Licensees" and
+"recipients" may be individuals or organizations.
+
+ To "modify" a work means to copy from or adapt all or part of the work
+in a fashion requiring copyright permission, other than the making of an
+exact copy. The resulting work is called a "modified version" of the
+earlier work or a work "based on" the earlier work.
+
+ A "covered work" means either the unmodified Program or a work based
+on the Program.
+
+ To "propagate" a work means to do anything with it that, without
+permission, would make you directly or secondarily liable for
+infringement under applicable copyright law, except executing it on a
+computer or modifying a private copy. Propagation includes copying,
+distribution (with or without modification), making available to the
+public, and in some countries other activities as well.
+
+ To "convey" a work means any kind of propagation that enables other
+parties to make or receive copies. Mere interaction with a user through
+a computer network, with no transfer of a copy, is not conveying.
+
+ An interactive user interface displays "Appropriate Legal Notices"
+to the extent that it includes a convenient and prominently visible
+feature that (1) displays an appropriate copyright notice, and (2)
+tells the user that there is no warranty for the work (except to the
+extent that warranties are provided), that licensees may convey the
+work under this License, and how to view a copy of this License. If
+the interface presents a list of user commands or options, such as a
+menu, a prominent item in the list meets this criterion.
+
+ 1. Source Code.
+
+ The "source code" for a work means the preferred form of the work
+for making modifications to it. "Object code" means any non-source
+form of a work.
+
+ A "Standard Interface" means an interface that either is an official
+standard defined by a recognized standards body, or, in the case of
+interfaces specified for a particular programming language, one that
+is widely used among developers working in that language.
+
+ The "System Libraries" of an executable work include anything, other
+than the work as a whole, that (a) is included in the normal form of
+packaging a Major Component, but which is not part of that Major
+Component, and (b) serves only to enable use of the work with that
+Major Component, or to implement a Standard Interface for which an
+implementation is available to the public in source code form. A
+"Major Component", in this context, means a major essential component
+(kernel, window system, and so on) of the specific operating system
+(if any) on which the executable work runs, or a compiler used to
+produce the work, or an object code interpreter used to run it.
+
+ The "Corresponding Source" for a work in object code form means all
+the source code needed to generate, install, and (for an executable
+work) run the object code and to modify the work, including scripts to
+control those activities. However, it does not include the work's
+System Libraries, or general-purpose tools or generally available free
+programs which are used unmodified in performing those activities but
+which are not part of the work. For example, Corresponding Source
+includes interface definition files associated with source files for
+the work, and the source code for shared libraries and dynamically
+linked subprograms that the work is specifically designed to require,
+such as by intimate data communication or control flow between those
+subprograms and other parts of the work.
+
+ The Corresponding Source need not include anything that users
+can regenerate automatically from other parts of the Corresponding
+Source.
+
+ The Corresponding Source for a work in source code form is that
+same work.
+
+ 2. Basic Permissions.
+
+ All rights granted under this License are granted for the term of
+copyright on the Program, and are irrevocable provided the stated
+conditions are met. This License explicitly affirms your unlimited
+permission to run the unmodified Program. The output from running a
+covered work is covered by this License only if the output, given its
+content, constitutes a covered work. This License acknowledges your
+rights of fair use or other equivalent, as provided by copyright law.
+
+ You may make, run and propagate covered works that you do not
+convey, without conditions so long as your license otherwise remains
+in force. You may convey covered works to others for the sole purpose
+of having them make modifications exclusively for you, or provide you
+with facilities for running those works, provided that you comply with
+the terms of this License in conveying all material for which you do
+not control copyright. Those thus making or running the covered works
+for you must do so exclusively on your behalf, under your direction
+and control, on terms that prohibit them from making any copies of
+your copyrighted material outside their relationship with you.
+
+ Conveying under any other circumstances is permitted solely under
+the conditions stated below. Sublicensing is not allowed; section 10
+makes it unnecessary.
+
+ 3. Protecting Users' Legal Rights From Anti-Circumvention Law.
+
+ No covered work shall be deemed part of an effective technological
+measure under any applicable law fulfilling obligations under article
+11 of the WIPO copyright treaty adopted on 20 December 1996, or
+similar laws prohibiting or restricting circumvention of such
+measures.
+
+ When you convey a covered work, you waive any legal power to forbid
+circumvention of technological measures to the extent such circumvention
+is effected by exercising rights under this License with respect to
+the covered work, and you disclaim any intention to limit operation or
+modification of the work as a means of enforcing, against the work's
+users, your or third parties' legal rights to forbid circumvention of
+technological measures.
+
+ 4. Conveying Verbatim Copies.
+
+ You may convey verbatim copies of the Program's source code as you
+receive it, in any medium, provided that you conspicuously and
+appropriately publish on each copy an appropriate copyright notice;
+keep intact all notices stating that this License and any
+non-permissive terms added in accord with section 7 apply to the code;
+keep intact all notices of the absence of any warranty; and give all
+recipients a copy of this License along with the Program.
+
+ You may charge any price or no price for each copy that you convey,
+and you may offer support or warranty protection for a fee.
+
+ 5. Conveying Modified Source Versions.
+
+ You may convey a work based on the Program, or the modifications to
+produce it from the Program, in the form of source code under the
+terms of section 4, provided that you also meet all of these conditions:
+
+ a) The work must carry prominent notices stating that you modified
+ it, and giving a relevant date.
+
+ b) The work must carry prominent notices stating that it is
+ released under this License and any conditions added under section
+ 7. This requirement modifies the requirement in section 4 to
+ "keep intact all notices".
+
+ c) You must license the entire work, as a whole, under this
+ License to anyone who comes into possession of a copy. This
+ License will therefore apply, along with any applicable section 7
+ additional terms, to the whole of the work, and all its parts,
+ regardless of how they are packaged. This License gives no
+ permission to license the work in any other way, but it does not
+ invalidate such permission if you have separately received it.
+
+ d) If the work has interactive user interfaces, each must display
+ Appropriate Legal Notices; however, if the Program has interactive
+ interfaces that do not display Appropriate Legal Notices, your
+ work need not make them do so.
+
+ A compilation of a covered work with other separate and independent
+works, which are not by their nature extensions of the covered work,
+and which are not combined with it such as to form a larger program,
+in or on a volume of a storage or distribution medium, is called an
+"aggregate" if the compilation and its resulting copyright are not
+used to limit the access or legal rights of the compilation's users
+beyond what the individual works permit. Inclusion of a covered work
+in an aggregate does not cause this License to apply to the other
+parts of the aggregate.
+
+ 6. Conveying Non-Source Forms.
+
+ You may convey a covered work in object code form under the terms
+of sections 4 and 5, provided that you also convey the
+machine-readable Corresponding Source under the terms of this License,
+in one of these ways:
+
+ a) Convey the object code in, or embodied in, a physical product
+ (including a physical distribution medium), accompanied by the
+ Corresponding Source fixed on a durable physical medium
+ customarily used for software interchange.
+
+ b) Convey the object code in, or embodied in, a physical product
+ (including a physical distribution medium), accompanied by a
+ written offer, valid for at least three years and valid for as
+ long as you offer spare parts or customer support for that product
+ model, to give anyone who possesses the object code either (1) a
+ copy of the Corresponding Source for all the software in the
+ product that is covered by this License, on a durable physical
+ medium customarily used for software interchange, for a price no
+ more than your reasonable cost of physically performing this
+ conveying of source, or (2) access to copy the
+ Corresponding Source from a network server at no charge.
+
+ c) Convey individual copies of the object code with a copy of the
+ written offer to provide the Corresponding Source. This
+ alternative is allowed only occasionally and noncommercially, and
+ only if you received the object code with such an offer, in accord
+ with subsection 6b.
+
+ d) Convey the object code by offering access from a designated
+ place (gratis or for a charge), and offer equivalent access to the
+ Corresponding Source in the same way through the same place at no
+ further charge. You need not require recipients to copy the
+ Corresponding Source along with the object code. If the place to
+ copy the object code is a network server, the Corresponding Source
+ may be on a different server (operated by you or a third party)
+ that supports equivalent copying facilities, provided you maintain
+ clear directions next to the object code saying where to find the
+ Corresponding Source. Regardless of what server hosts the
+ Corresponding Source, you remain obligated to ensure that it is
+ available for as long as needed to satisfy these requirements.
+
+ e) Convey the object code using peer-to-peer transmission, provided
+ you inform other peers where the object code and Corresponding
+ Source of the work are being offered to the general public at no
+ charge under subsection 6d.
+
+ A separable portion of the object code, whose source code is excluded
+from the Corresponding Source as a System Library, need not be
+included in conveying the object code work.
+
+ A "User Product" is either (1) a "consumer product", which means any
+tangible personal property which is normally used for personal, family,
+or household purposes, or (2) anything designed or sold for incorporation
+into a dwelling. In determining whether a product is a consumer product,
+doubtful cases shall be resolved in favor of coverage. For a particular
+product received by a particular user, "normally used" refers to a
+typical or common use of that class of product, regardless of the status
+of the particular user or of the way in which the particular user
+actually uses, or expects or is expected to use, the product. A product
+is a consumer product regardless of whether the product has substantial
+commercial, industrial or non-consumer uses, unless such uses represent
+the only significant mode of use of the product.
+
+ "Installation Information" for a User Product means any methods,
+procedures, authorization keys, or other information required to install
+and execute modified versions of a covered work in that User Product from
+a modified version of its Corresponding Source. The information must
+suffice to ensure that the continued functioning of the modified object
+code is in no case prevented or interfered with solely because
+modification has been made.
+
+ If you convey an object code work under this section in, or with, or
+specifically for use in, a User Product, and the conveying occurs as
+part of a transaction in which the right of possession and use of the
+User Product is transferred to the recipient in perpetuity or for a
+fixed term (regardless of how the transaction is characterized), the
+Corresponding Source conveyed under this section must be accompanied
+by the Installation Information. But this requirement does not apply
+if neither you nor any third party retains the ability to install
+modified object code on the User Product (for example, the work has
+been installed in ROM).
+
+ The requirement to provide Installation Information does not include a
+requirement to continue to provide support service, warranty, or updates
+for a work that has been modified or installed by the recipient, or for
+the User Product in which it has been modified or installed. Access to a
+network may be denied when the modification itself materially and
+adversely affects the operation of the network or violates the rules and
+protocols for communication across the network.
+
+ Corresponding Source conveyed, and Installation Information provided,
+in accord with this section must be in a format that is publicly
+documented (and with an implementation available to the public in
+source code form), and must require no special password or key for
+unpacking, reading or copying.
+
+ 7. Additional Terms.
+
+ "Additional permissions" are terms that supplement the terms of this
+License by making exceptions from one or more of its conditions.
+Additional permissions that are applicable to the entire Program shall
+be treated as though they were included in this License, to the extent
+that they are valid under applicable law. If additional permissions
+apply only to part of the Program, that part may be used separately
+under those permissions, but the entire Program remains governed by
+this License without regard to the additional permissions.
+
+ When you convey a copy of a covered work, you may at your option
+remove any additional permissions from that copy, or from any part of
+it. (Additional permissions may be written to require their own
+removal in certain cases when you modify the work.) You may place
+additional permissions on material, added by you to a covered work,
+for which you have or can give appropriate copyright permission.
+
+ Notwithstanding any other provision of this License, for material you
+add to a covered work, you may (if authorized by the copyright holders of
+that material) supplement the terms of this License with terms:
+
+ a) Disclaiming warranty or limiting liability differently from the
+ terms of sections 15 and 16 of this License; or
+
+ b) Requiring preservation of specified reasonable legal notices or
+ author attributions in that material or in the Appropriate Legal
+ Notices displayed by works containing it; or
+
+ c) Prohibiting misrepresentation of the origin of that material, or
+ requiring that modified versions of such material be marked in
+ reasonable ways as different from the original version; or
+
+ d) Limiting the use for publicity purposes of names of licensors or
+ authors of the material; or
+
+ e) Declining to grant rights under trademark law for use of some
+ trade names, trademarks, or service marks; or
+
+ f) Requiring indemnification of licensors and authors of that
+ material by anyone who conveys the material (or modified versions of
+ it) with contractual assumptions of liability to the recipient, for
+ any liability that these contractual assumptions directly impose on
+ those licensors and authors.
+
+ All other non-permissive additional terms are considered "further
+restrictions" within the meaning of section 10. If the Program as you
+received it, or any part of it, contains a notice stating that it is
+governed by this License along with a term that is a further
+restriction, you may remove that term. If a license document contains
+a further restriction but permits relicensing or conveying under this
+License, you may add to a covered work material governed by the terms
+of that license document, provided that the further restriction does
+not survive such relicensing or conveying.
+
+ If you add terms to a covered work in accord with this section, you
+must place, in the relevant source files, a statement of the
+additional terms that apply to those files, or a notice indicating
+where to find the applicable terms.
+
+ Additional terms, permissive or non-permissive, may be stated in the
+form of a separately written license, or stated as exceptions;
+the above requirements apply either way.
+
+ 8. Termination.
+
+ You may not propagate or modify a covered work except as expressly
+provided under this License. Any attempt otherwise to propagate or
+modify it is void, and will automatically terminate your rights under
+this License (including any patent licenses granted under the third
+paragraph of section 11).
+
+ However, if you cease all violation of this License, then your
+license from a particular copyright holder is reinstated (a)
+provisionally, unless and until the copyright holder explicitly and
+finally terminates your license, and (b) permanently, if the copyright
+holder fails to notify you of the violation by some reasonable means
+prior to 60 days after the cessation.
+
+ Moreover, your license from a particular copyright holder is
+reinstated permanently if the copyright holder notifies you of the
+violation by some reasonable means, this is the first time you have
+received notice of violation of this License (for any work) from that
+copyright holder, and you cure the violation prior to 30 days after
+your receipt of the notice.
+
+ Termination of your rights under this section does not terminate the
+licenses of parties who have received copies or rights from you under
+this License. If your rights have been terminated and not permanently
+reinstated, you do not qualify to receive new licenses for the same
+material under section 10.
+
+ 9. Acceptance Not Required for Having Copies.
+
+ You are not required to accept this License in order to receive or
+run a copy of the Program. Ancillary propagation of a covered work
+occurring solely as a consequence of using peer-to-peer transmission
+to receive a copy likewise does not require acceptance. However,
+nothing other than this License grants you permission to propagate or
+modify any covered work. These actions infringe copyright if you do
+not accept this License. Therefore, by modifying or propagating a
+covered work, you indicate your acceptance of this License to do so.
+
+ 10. Automatic Licensing of Downstream Recipients.
+
+ Each time you convey a covered work, the recipient automatically
+receives a license from the original licensors, to run, modify and
+propagate that work, subject to this License. You are not responsible
+for enforcing compliance by third parties with this License.
+
+ An "entity transaction" is a transaction transferring control of an
+organization, or substantially all assets of one, or subdividing an
+organization, or merging organizations. If propagation of a covered
+work results from an entity transaction, each party to that
+transaction who receives a copy of the work also receives whatever
+licenses to the work the party's predecessor in interest had or could
+give under the previous paragraph, plus a right to possession of the
+Corresponding Source of the work from the predecessor in interest, if
+the predecessor has it or can get it with reasonable efforts.
+
+ You may not impose any further restrictions on the exercise of the
+rights granted or affirmed under this License. For example, you may
+not impose a license fee, royalty, or other charge for exercise of
+rights granted under this License, and you may not initiate litigation
+(including a cross-claim or counterclaim in a lawsuit) alleging that
+any patent claim is infringed by making, using, selling, offering for
+sale, or importing the Program or any portion of it.
+
+ 11. Patents.
+
+ A "contributor" is a copyright holder who authorizes use under this
+License of the Program or a work on which the Program is based. The
+work thus licensed is called the contributor's "contributor version".
+
+ A contributor's "essential patent claims" are all patent claims
+owned or controlled by the contributor, whether already acquired or
+hereafter acquired, that would be infringed by some manner, permitted
+by this License, of making, using, or selling its contributor version,
+but do not include claims that would be infringed only as a
+consequence of further modification of the contributor version. For
+purposes of this definition, "control" includes the right to grant
+patent sublicenses in a manner consistent with the requirements of
+this License.
+
+ Each contributor grants you a non-exclusive, worldwide, royalty-free
+patent license under the contributor's essential patent claims, to
+make, use, sell, offer for sale, import and otherwise run, modify and
+propagate the contents of its contributor version.
+
+ In the following three paragraphs, a "patent license" is any express
+agreement or commitment, however denominated, not to enforce a patent
+(such as an express permission to practice a patent or covenant not to
+sue for patent infringement). To "grant" such a patent license to a
+party means to make such an agreement or commitment not to enforce a
+patent against the party.
+
+ If you convey a covered work, knowingly relying on a patent license,
+and the Corresponding Source of the work is not available for anyone
+to copy, free of charge and under the terms of this License, through a
+publicly available network server or other readily accessible means,
+then you must either (1) cause the Corresponding Source to be so
+available, or (2) arrange to deprive yourself of the benefit of the
+patent license for this particular work, or (3) arrange, in a manner
+consistent with the requirements of this License, to extend the patent
+license to downstream recipients. "Knowingly relying" means you have
+actual knowledge that, but for the patent license, your conveying the
+covered work in a country, or your recipient's use of the covered work
+in a country, would infringe one or more identifiable patents in that
+country that you have reason to believe are valid.
+
+ If, pursuant to or in connection with a single transaction or
+arrangement, you convey, or propagate by procuring conveyance of, a
+covered work, and grant a patent license to some of the parties
+receiving the covered work authorizing them to use, propagate, modify
+or convey a specific copy of the covered work, then the patent license
+you grant is automatically extended to all recipients of the covered
+work and works based on it.
+
+ A patent license is "discriminatory" if it does not include within
+the scope of its coverage, prohibits the exercise of, or is
+conditioned on the non-exercise of one or more of the rights that are
+specifically granted under this License. You may not convey a covered
+work if you are a party to an arrangement with a third party that is
+in the business of distributing software, under which you make payment
+to the third party based on the extent of your activity of conveying
+the work, and under which the third party grants, to any of the
+parties who would receive the covered work from you, a discriminatory
+patent license (a) in connection with copies of the covered work
+conveyed by you (or copies made from those copies), or (b) primarily
+for and in connection with specific products or compilations that
+contain the covered work, unless you entered into that arrangement,
+or that patent license was granted, prior to 28 March 2007.
+
+ Nothing in this License shall be construed as excluding or limiting
+any implied license or other defenses to infringement that may
+otherwise be available to you under applicable patent law.
+
+ 12. No Surrender of Others' Freedom.
+
+ If conditions are imposed on you (whether by court order, agreement or
+otherwise) that contradict the conditions of this License, they do not
+excuse you from the conditions of this License. If you cannot convey a
+covered work so as to satisfy simultaneously your obligations under this
+License and any other pertinent obligations, then as a consequence you may
+not convey it at all. For example, if you agree to terms that obligate you
+to collect a royalty for further conveying from those to whom you convey
+the Program, the only way you could satisfy both those terms and this
+License would be to refrain entirely from conveying the Program.
+
+ 13. Use with the GNU Affero General Public License.
+
+ Notwithstanding any other provision of this License, you have
+permission to link or combine any covered work with a work licensed
+under version 3 of the GNU Affero General Public License into a single
+combined work, and to convey the resulting work. The terms of this
+License will continue to apply to the part which is the covered work,
+but the special requirements of the GNU Affero General Public License,
+section 13, concerning interaction through a network will apply to the
+combination as such.
+
+ 14. Revised Versions of this License.
+
+ The Free Software Foundation may publish revised and/or new versions of
+the GNU General Public License from time to time. Such new versions will
+be similar in spirit to the present version, but may differ in detail to
+address new problems or concerns.
+
+ Each version is given a distinguishing version number. If the
+Program specifies that a certain numbered version of the GNU General
+Public License "or any later version" applies to it, you have the
+option of following the terms and conditions either of that numbered
+version or of any later version published by the Free Software
+Foundation. If the Program does not specify a version number of the
+GNU General Public License, you may choose any version ever published
+by the Free Software Foundation.
+
+ If the Program specifies that a proxy can decide which future
+versions of the GNU General Public License can be used, that proxy's
+public statement of acceptance of a version permanently authorizes you
+to choose that version for the Program.
+
+ Later license versions may give you additional or different
+permissions. However, no additional obligations are imposed on any
+author or copyright holder as a result of your choosing to follow a
+later version.
+
+ 15. Disclaimer of Warranty.
+
+ THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY
+APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT
+HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY
+OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO,
+THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM
+IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF
+ALL NECESSARY SERVICING, REPAIR OR CORRECTION.
+
+ 16. Limitation of Liability.
+
+ IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS
+THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY
+GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE
+USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF
+DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD
+PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS),
+EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF
+SUCH DAMAGES.
+
+ 17. Interpretation of Sections 15 and 16.
+
+ If the disclaimer of warranty and limitation of liability provided
+above cannot be given local legal effect according to their terms,
+reviewing courts shall apply local law that most closely approximates
+an absolute waiver of all civil liability in connection with the
+Program, unless a warranty or assumption of liability accompanies a
+copy of the Program in return for a fee.
+
+ END OF TERMS AND CONDITIONS
+
+ How to Apply These Terms to Your New Programs
+
+ If you develop a new program, and you want it to be of the greatest
+possible use to the public, the best way to achieve this is to make it
+free software which everyone can redistribute and change under these terms.
+
+ To do so, attach the following notices to the program. It is safest
+to attach them to the start of each source file to most effectively
+state the exclusion of warranty; and each file should have at least
+the "copyright" line and a pointer to where the full notice is found.
+
+ <one line to give the program's name and a brief idea of what it does.>
+ Copyright (C) <year> <name of author>
+
+ This program is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+Also add information on how to contact you by electronic and paper mail.
+
+ If the program does terminal interaction, make it output a short
+notice like this when it starts in an interactive mode:
+
+ <program> Copyright (C) <year> <name of author>
+ This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
+ This is free software, and you are welcome to redistribute it
+ under certain conditions; type `show c' for details.
+
+The hypothetical commands `show w' and `show c' should show the appropriate
+parts of the General Public License. Of course, your program's commands
+might be different; for a GUI interface, you would use an "about box".
+
+ You should also get your employer (if you work as a programmer) or school,
+if any, to sign a "copyright disclaimer" for the program, if necessary.
+For more information on this, and how to apply and follow the GNU GPL, see
+<http://www.gnu.org/licenses/>.
+
+ The GNU General Public License does not permit incorporating your program
+into proprietary programs. If your program is a subroutine library, you
+may consider it more useful to permit linking proprietary applications with
+the library. If this is what you want to do, use the GNU Lesser General
+Public License instead of this License. But first, please read
+<http://www.gnu.org/philosophy/why-not-lgpl.html>.
Added: seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/README.md
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/README.md (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/README.md 2013-09-26 00:57:45 UTC (rev 22850)
@@ -0,0 +1,4 @@
+GEOCUBIT--experimental
+======================
+
+experimental version of GEOCUBIT
\ No newline at end of file
Modified: seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/doc/geocubit_manual.pdf
===================================================================
(Binary files differ)
Deleted: seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/doc/geocubit_manual.pdfsync
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/doc/geocubit_manual.pdfsync 2013-09-26 00:49:25 UTC (rev 22849)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/doc/geocubit_manual.pdfsync 2013-09-26 00:57:45 UTC (rev 22850)
@@ -1,1788 +0,0 @@
-geocubit_manual
-version 1
-(nameref.sty
-(refcount.sty
-)
-(gettitlestring.sty
-(gettitlestring.cfg
-)
-)
-)
-(.geocubit_manual.out
-)
-(.geocubit_manual.out
-)
-l 0 41
-l 1 41
-l 2 41
-l 3 41
-l 4 41
-l 5 41
-l 6 41
-l 7 41
-l 8 41
-l 9 41
-l 10 41
-l 11 41
-l 12 41
-l 13 41
-l 14 41
-l 15 42
-l 16 42
-l 17 42
-l 18 42
-l 19 42
-l 20 42
-l 21 2
-l 22 2
-l 23 2
-l 24 3
-l 25 3
-l 26 3
-l 27 4
-l 28 4
-l 29 4
-l 30 4
-l 31 4
-l 32 4
-l 33 4
-l 34 5
-l 35 5
-l 36 5
-l 37 5
-l 38 5
-l 39 5
-l 40 5
-l 41 6
-l 42 6
-l 43 6
-l 44 7
-l 45 7
-l 46 7
-l 47 7
-l 48 7
-l 49 7
-l 50 7
-l 51 8
-l 52 8
-l 53 8
-l 54 8
-l 55 8
-l 56 8
-l 57 8
-l 58 9
-l 59 9
-l 60 9
-l 61 9
-l 62 9
-l 63 9
-l 64 9
-l 65 10
-l 66 10
-l 67 10
-l 68 10
-l 69 10
-l 70 10
-l 71 10
-l 72 11
-l 73 11
-l 74 11
-l 75 12
-l 76 12
-l 77 12
-l 78 13
-l 79 13
-l 80 13
-l 81 13
-l 82 13
-l 83 13
-l 84 13
-l 85 14
-l 86 14
-l 87 14
-l 88 14
-l 89 14
-l 90 14
-l 91 14
-l 92 15
-l 93 15
-l 94 15
-l 95 15
-l 96 15
-l 97 15
-l 98 15
-l 99 16
-l 100 16
-l 101 16
-l 102 16
-l 103 16
-l 104 16
-l 105 16
-l 106 17
-l 107 17
-l 108 17
-l 109 17
-l 110 17
-l 111 17
-l 112 17
-l 113 18
-l 114 18
-l 115 18
-l 116 18
-l 117 18
-l 118 18
-l 119 18
-l 120 19
-l 121 19
-l 122 19
-l 123 19
-l 124 19
-l 125 19
-l 126 19
-l 127 20
-l 128 20
-l 129 20
-l 130 20
-l 131 20
-l 132 20
-l 133 20
-l 134 21
-l 135 21
-l 136 21
-l 137 22
-l 138 22
-l 139 22
-l 140 43
-l 141 43
-l 142 43
-l 143 43
-l 144 43
-l 145 43
-l 146 43
-l 147 43
-l 148 43
-l 149 43
-l 150 43
-l 151 43
-l 152 43
-l 153 43
-s 1
-p 153 4736286 50644704
-p 149 4736286 50644704
-p 150 4736286 50644704
-p 151 4736286 49530837
-p 141 4736286 50644704
-p 143 4736286 49465056
-p 144 4736286 49465056
-p 145 4736286 48678624
-p 140 4736286 47040224
-p 5 12334523 43450440
-p 2 12334523 43450440
-p 3 11975713 43450440
-p 1 11975713 43450440
-p 7 15755887 41462614
-p 8 15755887 41462614
-p 9 15755887 41462614
-p 11 15755887 42104864
-p 10 15755887 41462614
-p 12 16149103 41462614
-p 14 16450632 39857776
-p 17 4736286 36810553
-p 16 4736286 36810553
-p 18 4736286 36810553
-p 19 4736286 36810553
-p 20 4736286 37990201
-p 21 5806608 35205715
-p 22 4736286 35205715
-p 23 32956737 35205715
-p 24 5806608 33600877
-p 25 4736286 33600877
-p 26 32956737 33600877
-p 27 7447770 32709587
-p 29 5806608 32709587
-p 33 33121485 32709587
-p 34 7447770 31818297
-p 36 5806608 31818297
-p 40 33121485 31818297
-p 41 5806608 30213459
-p 42 4736286 30213459
-p 43 32956737 30213459
-p 44 7447770 29322169
-p 46 5806608 29322169
-p 50 33121485 29322169
-p 51 7447770 28430879
-p 53 5806608 28430879
-p 57 33121485 28430879
-p 58 7447770 27539589
-p 60 5806608 27539589
-p 64 33121485 27539589
-p 65 7447770 26648299
-p 67 5806608 26648299
-p 71 33121485 26648299
-p 72 5806608 25043461
-p 73 4736286 25043461
-p 74 32956737 25043461
-p 75 5806608 23438623
-p 76 4736286 23438623
-p 77 32956737 23438623
-p 78 7447770 22547333
-p 80 5806608 22547333
-p 84 33121485 22547333
-p 85 7447770 21656043
-p 87 5806608 21656043
-p 91 33121485 21656043
-p 92 7447770 20764753
-p 94 5806608 20764753
-p 98 33121485 20764753
-p 99 7447770 19873463
-p 101 5806608 19873463
-p 105 33121485 19873463
-p 106 7447770 18982173
-p 108 5806608 18982173
-p 112 33121485 18982173
-p 113 7447770 18090883
-p 115 5806608 18090883
-p 119 33121485 18090883
-p 120 7447770 17199593
-p 122 5806608 17199593
-p 126 33121485 17199593
-p 127 7447770 16308303
-p 129 5806608 16308303
-p 133 33121485 16308303
-p 134 5806608 14703465
-p 135 4736286 14703465
-p 136 32956737 14703465
-p 137 5806608 13098627
-p 138 4736286 13098627
-p 139 32956737 11316047
-p 147 4736286 6027778
-p 148 4736286 6027778
-l 154 45
-l 155 45
-l 156 45
-l 157 47
-l 158 52
-l 159 52
-(omscmr.fd
-)
-l 160 52
-l 161 52
-l 162 52
-l 163 53
-l 164 53
-l 165 53
-l 166 53
-l 167 53
-l 168 54
-l 169 54
-l 170 54
-l 171 54
-l 172 54
-l 173 56
-l 174 56
-l 175 56
-l 176 56
-l 177 56
-l 178 60
-l 179 64
-l 180 64
-l 181 64
-l 182 64
-l 183 64
-l 184 64
-l 185 64
-l 186 64
-l 187 64
-l 188 64
-l 189 64
-l 190 64
-l 191 64
-s 2
-p 191 4736286 6027778
-p 188 4736286 50644704
-p 189 4736286 50644704
-p 190 4736286 49530837
-p 180 4736286 50644704
-p 182 4736286 49465056
-p 183 4736286 49465056
-p 184 4736286 48678624
-p 179 4736286 47040224
-p 156 4736286 46319328
-p 155 4736286 46319328
-p 157 4736286 44717572
-p 162 6530338 39474690
-p 160 6530338 39474690
-p 161 4736286 39474690
-p 158 4736286 39474690
-p 159 6171528 39474690
-p 167 6530338 37993576
-p 165 6530338 37993576
-p 166 4736286 37993576
-p 163 4736286 37993576
-p 164 6171528 37993576
-p 172 6530338 36512462
-p 170 6530338 36512462
-p 171 4736286 36512462
-p 168 4736286 36512462
-p 169 6171528 36512462
-p 177 6530338 34140058
-p 175 6530338 34140058
-p 176 4736286 34140058
-p 173 4736286 34140058
-p 174 6171528 34140058
-p 178 4736286 30679756
-p 186 4736286 6027778
-p 187 4736286 6027778
-l 192 65
-l 193 65
-l 194 65
-l 195 68
-l 196 68
-l 197 68
-l 198 70
-l 199 73
-l 200 73
-l 201 73
-l 202 73
-l 203 73
-l 204 74
-l 205 74
-l 206 74
-l 207 74
-l 208 74
-l 209 75
-l 210 75
-l 211 75
-l 212 75
-l 213 75
-l 214 78
-l 215 80
-l 216 80
-l 217 80
-l 218 80
-l 219 80
-l 220 83
-l 221 83
-l 222 83
-l 223 85
-(t1cmtt.fd
-)
-l 224 86
-l 225 86
-l 226 86
-l 227 87
-l 228 90
-l 229 92
-l 230 92
-l 231 92
-l 232 93
-l 233 94
-l 234 95
-l 235 96
-l 236 97
-l 237 101
-l 238 101
-l 239 101
-l 240 101
-l 241 101
-l 242 101
-l 243 101
-l 244 101
-l 245 101
-l 246 101
-l 247 101
-l 248 101
-l 249 101
-s 3
-p 249 4736286 6027778
-p 246 4736286 50644704
-p 247 4736286 50644704
-p 248 4736286 49530837
-p 238 4736286 50644704
-p 240 4736286 49465056
-p 241 4736286 49465056
-p 242 4736286 48678624
-p 237 4736286 47040224
-p 194 4736286 46319328
-p 193 4736286 46319328
-p 197 4736286 44691358
-p 196 4736286 44691358
-p 198 4736286 43336721
-p 203 6530338 41658999
-p 201 6530338 41658999
-p 202 4736286 41658999
-p 199 4736286 41658999
-p 200 6171528 41658999
-p 208 6530338 40177885
-p 206 6530338 40177885
-p 207 4736286 40177885
-p 204 4736286 40177885
-p 205 6171528 40177885
-p 213 6530338 38696771
-p 211 6530338 38696771
-p 212 4736286 38696771
-p 209 4736286 38696771
-p 210 6171528 38696771
-p 214 4736286 37019049
-p 219 6530338 35537935
-p 217 6530338 35537935
-p 218 4736286 35537935
-p 215 4736286 35537935
-p 216 6171528 35537935
-p 222 4736286 33616513
-p 221 4736286 33616513
-p 223 4736286 32261876
-p 227 6530338 30584154
-p 225 6530338 30584154
-p 226 4736286 30584154
-p 224 4736286 30584154
-p 228 4736286 28906432
-p 232 6530338 27228710
-p 230 6530338 27228710
-p 231 4736286 27228710
-p 229 4736286 27228710
-p 233 6530338 26337420
-p 234 6530338 25446130
-p 235 6530338 24554840
-p 236 6530338 23663550
-p 244 4736286 6027778
-p 245 4736286 6027778
-l 250 102
-l 251 102
-l 252 102
-l 253 104
-l 254 104
-l 255 104
-l 256 106
-l 257 106
-l 258 106
-l 259 106
-l 260 106
-l 261 107
-l 262 107
-l 263 107
-l 264 108
-l 265 111
-l 266 111
-l 267 111
-l 268 111
-l 269 111
-l 270 112
-l 271 112
-l 272 112
-l 273 113
-l 274 117
-l 275 117
-l 276 117
-l 277 118
-l 278 121
-l 279 121
-l 280 121
-l 281 121
-l 282 121
-l 283 122
-l 284 122
-l 285 122
-l 286 123
-l 287 128
-l 288 128
-l 289 128
-l 290 128
-l 291 128
-l 292 129
-l 293 129
-l 294 129
-l 295 130
-l 296 134
-l 297 136
-l 298 136
-l 299 136
-l 300 136
-l 301 136
-l 302 137
-l 303 137
-l 304 137
-l 305 138
-l 306 143
-l 307 143
-l 308 143
-l 309 143
-l 310 143
-l 311 144
-l 312 144
-l 313 144
-l 314 145
-l 315 149
-l 316 149
-l 317 149
-l 318 149
-l 319 149
-l 320 149
-l 321 149
-l 322 151
-l 323 151
-l 324 151
-l 325 151
-l 326 151
-l 327 151
-l 328 151
-l 329 151
-l 330 151
-l 331 151
-l 332 151
-l 333 151
-l 334 151
-l 335 151
-l 336 152
-l 337 152
-l 338 153
-l 339 153
-l 340 153
-l 341 153
-l 342 153
-l 343 153
-l 344 153
-l 345 153
-l 346 153
-l 347 153
-l 348 153
-l 349 153
-l 350 153
-l 351 154
-l 352 154
-l 353 154
-l 354 154
-l 355 154
-l 356 154
-l 357 154
-l 358 157
-l 359 157
-l 360 157
-l 361 159
-l 362 162
-l 363 162
-l 364 162
-l 365 164
-l 366 164
-l 367 164
-l 368 164
-l 369 164
-l 370 165
-l 371 165
-l 372 165
-l 373 166
-l 374 171
-l 375 171
-l 376 171
-l 377 171
-l 378 171
-l 379 172
-l 380 172
-l 381 172
-l 382 173
-l 383 175
-l 384 175
-l 385 179
-l 386 179
-l 387 179
-l 388 179
-l 389 179
-l 390 180
-l 391 180
-l 392 180
-l 393 181
-l 394 182
-l 395 182
-l 396 182
-l 397 182
-l 398 182
-l 399 182
-l 400 182
-l 401 182
-l 402 182
-l 403 182
-l 404 182
-l 405 182
-l 406 182
-s 4
-p 406 4736286 6027778
-p 403 4736286 50644704
-p 404 4736286 50644704
-p 405 4736286 49530837
-p 395 4736286 50644704
-p 397 4736286 49465056
-p 398 4736286 49465056
-p 399 4736286 48678624
-p 394 4736286 47040224
-p 360 4736286 46319328
-p 358 4736286 46319328
-p 252 4736286 45598432
-p 251 4736286 45598432
-p 255 4736286 43970462
-p 254 4736286 43970462
-p 260 6530338 42615825
-p 258 6530338 42615825
-p 259 4736286 42615825
-p 256 4736286 42615825
-p 257 6171528 42615825
-p 264 8109102 40938103
-p 262 8109102 40938103
-p 263 6530338 40938103
-p 261 6530338 40938103
-p 269 6530338 39260381
-p 267 6530338 39260381
-p 268 4736286 39260381
-p 265 4736286 39260381
-p 266 6171528 39260381
-p 273 8109102 37582659
-p 271 8109102 37582659
-p 272 6530338 37582659
-p 270 6530338 37582659
-p 276 4736286 35661237
-p 275 4736286 35661237
-p 277 4736286 34306600
-p 282 6530338 32628878
-p 280 6530338 32628878
-p 281 4736286 32628878
-p 278 4736286 32628878
-p 279 6171528 32628878
-p 286 8109102 30951156
-p 284 8109102 30951156
-p 285 6530338 30951156
-p 283 6530338 30951156
-p 291 6530338 29273434
-p 289 6530338 29273434
-p 290 4736286 29273434
-p 287 4736286 29273434
-p 288 6171528 29273434
-p 295 8109102 27595712
-p 293 8109102 27595712
-p 294 6530338 27595712
-p 292 6530338 27595712
-p 296 4736286 25917990
-p 301 6530338 24436876
-p 299 6530338 24436876
-p 300 4736286 24436876
-p 297 4736286 24436876
-p 298 6171528 24436876
-p 305 8109102 22759154
-p 303 8109102 22759154
-p 304 6530338 22759154
-p 302 6530338 22759154
-p 310 6530338 21081432
-p 308 6530338 21081432
-p 309 4736286 21081432
-p 306 4736286 21081432
-p 307 6171528 21081432
-p 314 8109102 19403710
-p 312 8109102 19403710
-p 313 6530338 19403710
-p 311 6530338 19403710
-p 361 6530338 17725988
-p 364 4736286 14913276
-p 363 4736286 14913276
-p 369 6530338 13558639
-p 367 6530338 13558639
-p 368 4736286 13558639
-p 365 4736286 13558639
-p 366 6171528 13558639
-p 373 8109102 11880917
-p 371 8109102 11880917
-p 372 6530338 11880917
-p 370 6530338 11880917
-p 378 6530338 10203195
-p 376 6530338 10203195
-p 377 4736286 10203195
-p 374 4736286 10203195
-p 375 6171528 10203195
-p 382 8109102 8525473
-p 380 8109102 8525473
-p 381 6530338 8525473
-p 379 6530338 8525473
-p 401 4736286 6027778
-p 402 4736286 6027778
-l 407 185
-l 408 185
-l 409 185
-l 410 187
-l 411 190
-l 412 190
-l 413 190
-l 414 190
-l 415 190
-l 416 191
-l 417 191
-l 418 191
-l 419 192
-l 420 194
-l 421 195
-l 422 197
-l 423 197
-l 424 197
-l 425 197
-l 426 197
-l 427 198
-l 428 198
-l 429 198
-l 430 199
-l 431 201
-l 432 201
-l 433 205
-l 434 205
-l 435 205
-l 436 205
-l 437 205
-l 438 205
-l 439 205
-l 440 206
-l 441 206
-l 442 206
-l 443 206
-l 444 206
-l 445 206
-l 446 206
-l 447 206
-l 448 207
-l 449 207
-l 450 207
-l 451 207
-l 452 207
-l 453 207
-l 454 207
-l 455 208
-l 456 208
-l 457 208
-l 458 208
-l 459 208
-l 460 211
-l 461 211
-l 462 211
-l 463 211
-l 464 211
-l 465 211
-l 466 211
-l 467 211
-l 468 211
-l 469 211
-l 470 211
-l 471 211
-l 472 211
-l 473 211
-l 474 211
-l 475 211
-l 476 211
-l 477 211
-l 478 211
-s 5
-p 478 4736286 6027778
-p 475 4736286 50644704
-p 476 4736286 50644704
-p 477 4736286 49530837
-p 467 4736286 50644704
-p 469 4736286 49465056
-p 470 4736286 49465056
-p 471 4736286 48678624
-p 466 4736286 47040224
-p 465 4736286 47040224
-p 464 4736286 47040224
-p 315 4736286 47040224
-p 316 4736286 47040224
-p 322 8574223 38128105
-p 320 8574223 38128105
-p 321 8215413 38128105
-p 319 8215413 38128105
-p 331 8574223 45791660
-p 332 8574223 45791660
-p 323 8574223 38128105
-p 324 8574223 38128105
-p 326 8574223 38128105
-p 327 8574223 38128105
-p 329 8574223 38128105
-p 330 8574223 38128105
-p 328 8574223 38128105
-p 335 8574223 37236815
-p 333 12612061 37236815
-p 336 17659778 38128105
-p 337 17659778 39019395
-p 346 17659778 47040224
-p 347 17659778 47040224
-p 338 17659778 38128105
-p 339 17659778 38128105
-p 341 17659778 38128105
-p 342 17659778 38128105
-p 344 17659778 38128105
-p 345 17659778 38128105
-p 343 17659778 38128105
-p 350 17659778 37236815
-p 348 23502889 37236815
-p 355 4736286 35502730
-p 356 8022359 35502730
-p 357 8022359 36394020
-p 462 4736286 31720892
-p 460 4736286 31720892
-p 384 6530338 30999996
-p 389 6530338 29482571
-p 387 6530338 29482571
-p 388 4736286 29482571
-p 385 4736286 29482571
-p 386 6171528 29482571
-p 393 8109102 27884998
-p 391 8109102 27884998
-p 392 6530338 27884998
-p 390 6530338 27884998
-p 409 4736286 25097470
-p 408 4736286 25097470
-p 410 4736286 23742833
-p 415 6530338 21504010
-p 413 6530338 21504010
-p 414 4736286 21504010
-p 411 4736286 21504010
-p 412 6171528 21504010
-p 419 8109102 19015147
-p 417 8109102 19015147
-p 418 6530338 19015147
-p 416 6530338 19015147
-p 420 6530338 16526284
-p 421 6530338 14743704
-p 426 6530338 11533442
-p 424 6530338 11533442
-p 425 4736286 11533442
-p 422 4736286 11533442
-p 423 6171528 11533442
-p 430 8109102 9044579
-p 428 8109102 9044579
-p 429 6530338 9044579
-p 427 6530338 9044579
-p 473 4736286 6027778
-p 474 4736286 6027778
-l 479 214
-l 480 214
-l 481 214
-l 482 214
-l 483 214
-l 484 215
-l 485 215
-l 486 215
-l 487 216
-l 488 218
-l 489 218
-l 490 221
-l 491 221
-l 492 221
-l 493 221
-l 494 221
-l 495 221
-l 496 221
-l 497 221
-l 498 221
-l 499 221
-l 500 221
-l 501 221
-l 502 221
-l 503 221
-l 504 221
-l 505 221
-s 6
-p 505 4736286 6027778
-p 502 4736286 50644704
-p 503 4736286 50644704
-p 504 4736286 49530837
-p 494 4736286 50644704
-p 496 4736286 49465056
-p 497 4736286 49465056
-p 498 4736286 48678624
-p 493 4736286 47040224
-p 492 4736286 47040224
-p 491 4736286 47040224
-p 433 4736286 47040224
-p 434 4736286 47040224
-p 440 8695633 32629812
-p 438 8695633 32629812
-p 439 8336823 32629812
-p 437 8336823 32629812
-p 441 8695633 32629812
-p 443 8695633 32629812
-p 444 8695633 32629812
-p 446 8695633 32629812
-p 447 8695633 32629812
-p 445 8695633 32629812
-p 448 18528536 32629812
-p 450 18528536 32629812
-p 451 18528536 32629812
-p 453 18528536 32629812
-p 454 18528536 32629812
-p 452 18528536 32629812
-p 459 4736286 31083162
-p 456 17920895 31083162
-p 457 21042878 31083162
-p 458 21042878 31974452
-p 432 6530338 28912044
-p 483 6530338 26644498
-p 481 6530338 26644498
-p 482 4736286 26644498
-p 479 4736286 26644498
-p 480 6171528 26644498
-p 487 8109102 24966776
-p 485 8109102 24966776
-p 486 6530338 24966776
-p 484 6530338 24966776
-p 489 6530338 21676868
-p 500 4736286 6027778
-p 501 4736286 6027778
-l 506 222
-l 507 222
-l 508 222
-l 509 223
-l 510 225
-l 511 225
-l 512 225
-l 513 226
-l 514 226
-l 515 226
-l 516 229
-l 517 229
-l 518 229
-l 519 230
-l 520 230
-l 521 230
-l 522 231
-l 523 232
-l 524 235
-l 525 235
-l 526 235
-l 527 236
-l 528 236
-l 529 236
-l 530 237
-l 531 237
-l 532 237
-l 533 238
-l 534 238
-l 535 238
-l 536 239
-l 537 239
-l 538 239
-l 539 240
-l 540 240
-l 541 240
-l 542 242
-l 543 245
-l 544 245
-l 545 245
-l 546 246
-l 547 248
-l 548 248
-l 549 248
-l 550 250
-l 551 253
-l 552 255
-l 553 255
-l 554 255
-l 555 256
-l 556 257
-l 557 257
-l 558 257
-l 559 258
-l 560 258
-l 561 258
-l 562 259
-l 563 259
-l 564 259
-l 565 260
-l 566 260
-l 567 260
-l 568 261
-l 569 261
-l 570 261
-l 571 262
-l 572 262
-l 573 262
-l 574 263
-l 575 263
-l 576 263
-l 577 264
-l 578 264
-l 579 264
-l 580 265
-l 581 265
-l 582 265
-l 583 267
-l 584 267
-l 585 267
-l 586 268
-l 587 269
-l 588 269
-l 589 269
-l 590 273
-l 591 273
-l 592 273
-l 593 274
-l 594 274
-l 595 274
-l 596 274
-l 597 274
-l 598 274
-l 599 274
-l 600 274
-l 601 274
-l 602 274
-l 603 274
-l 604 274
-l 605 274
-s 7
-p 605 4736286 6027778
-p 602 4736286 50644704
-p 603 4736286 50644704
-p 604 4736286 49530837
-p 594 4736286 50644704
-p 596 4736286 49465056
-p 597 4736286 49465056
-p 598 4736286 48678624
-p 593 4736286 47040224
-p 508 4736286 46319328
-p 507 4736286 46319328
-p 509 4736286 44717572
-p 513 6530338 43039850
-p 511 6530338 43039850
-p 512 4736286 43039850
-p 510 4736286 43039850
-p 514 12463570 43039850
-p 519 6530338 41362128
-p 517 6530338 41362128
-p 518 4736286 41362128
-p 516 4736286 41362128
-p 520 18767629 41362128
-p 522 6530338 40470838
-p 523 6530338 39579548
-p 527 6530338 37901826
-p 525 6530338 37901826
-p 526 4736286 37901826
-p 524 4736286 37901826
-p 528 13576051 37901826
-p 530 6530338 37010536
-p 531 13946878 37010536
-p 536 6530338 35332814
-p 534 6530338 35332814
-p 535 4736286 35332814
-p 533 4736286 35332814
-p 537 13946878 35332814
-p 539 6530338 34441524
-p 540 19509283 34441524
-p 542 4736286 31872512
-p 545 4736286 29611721
-p 544 4736286 29611721
-p 546 4736286 28009965
-p 549 4736286 25197253
-p 548 4736286 25197253
-p 550 4736286 23842616
-p 551 4736286 22060036
-p 555 6530338 19491024
-p 553 6530338 19491024
-p 554 4736286 19491024
-p 552 4736286 19491024
-p 556 6530338 18599734
-p 557 12092743 18599734
-p 559 6530338 17708444
-p 560 11721916 17708444
-p 562 6530338 16817154
-p 563 12092743 16817154
-p 565 6530338 15925864
-p 566 11721916 15925864
-p 568 6530338 15034574
-p 569 12463570 15034574
-p 571 6530338 14143284
-p 572 13205224 14143284
-p 574 6530338 13251994
-p 575 16542667 13251994
-p 577 6530338 12360704
-p 578 16542667 12360704
-p 580 6530338 11469414
-p 581 16171840 11469414
-p 600 4736286 6027778
-p 601 4736286 6027778
-l 606 275
-l 607 275
-l 608 275
-l 609 276
-l 610 276
-l 611 276
-l 612 276
-l 613 277
-l 614 278
-l 615 278
-l 616 278
-l 617 280
-l 618 280
-l 619 280
-l 620 281
-l 621 282
-l 622 282
-l 623 284
-l 624 284
-l 625 284
-l 626 285
-l 627 286
-l 628 286
-l 629 286
-l 630 287
-l 631 287
-l 632 287
-l 633 288
-l 634 288
-l 635 288
-l 636 289
-l 637 289
-l 638 289
-l 639 290
-l 640 291
-l 641 291
-l 642 291
-l 643 294
-l 644 294
-l 645 294
-l 646 295
-l 647 295
-l 648 295
-l 649 296
-l 650 296
-l 651 297
-l 652 297
-l 653 297
-l 654 298
-l 655 298
-l 656 298
-l 657 299
-l 658 299
-l 659 299
-l 660 300
-l 661 303
-l 662 303
-l 663 303
-l 664 303
-l 665 303
-l 666 303
-l 667 303
-l 668 305
-l 669 305
-l 670 305
-l 671 305
-l 672 305
-l 673 305
-l 674 305
-l 675 305
-l 676 305
-l 677 305
-l 678 305
-l 679 305
-l 680 305
-l 681 305
-l 682 306
-l 683 306
-l 684 307
-l 685 307
-l 686 307
-l 687 307
-l 688 307
-l 689 307
-l 690 307
-l 691 307
-l 692 307
-l 693 307
-l 694 307
-l 695 307
-l 696 307
-l 697 308
-l 698 308
-l 699 309
-l 700 309
-l 701 309
-l 702 309
-l 703 309
-l 704 309
-l 705 309
-l 706 309
-l 707 309
-l 708 309
-l 709 309
-l 710 309
-l 711 309
-l 712 310
-l 713 310
-l 714 310
-l 715 310
-l 716 310
-l 717 310
-l 718 310
-l 719 313
-l 720 313
-l 721 313
-l 722 316
-l 723 316
-l 724 316
-l 725 317
-l 726 317
-l 727 317
-l 728 317
-l 729 318
-l 730 319
-l 731 321
-l 732 321
-l 733 321
-l 734 321
-l 735 321
-l 736 321
-l 737 321
-l 738 322
-l 739 322
-l 740 322
-l 741 322
-l 742 322
-l 743 322
-l 744 322
-l 745 322
-l 746 323
-l 747 323
-l 748 323
-l 749 323
-l 750 323
-l 751 323
-l 752 323
-l 753 324
-l 754 324
-l 755 324
-l 756 324
-l 757 324
-l 758 328
-l 759 328
-l 760 328
-l 761 332
-l 762 332
-l 763 332
-l 764 333
-l 765 334
-l 766 334
-l 767 334
-l 768 334
-l 769 334
-l 770 334
-l 771 334
-l 772 334
-l 773 334
-l 774 334
-l 775 334
-l 776 334
-l 777 334
-s 8
-p 777 4736286 6027778
-p 774 4736286 50644704
-p 775 4736286 50644704
-p 776 4736286 49530837
-p 766 4736286 50644704
-p 768 4736286 49465056
-p 769 4736286 49465056
-p 770 4736286 48678624
-p 765 4736286 47040224
-p 760 4736286 46319328
-p 758 4736286 46319328
-p 721 4736286 45598432
-p 719 4736286 45598432
-p 586 6530338 44877536
-p 584 6530338 44877536
-p 585 4736286 44877536
-p 583 4736286 44877536
-p 587 6530338 43986246
-p 588 10238608 43986246
-p 591 4736286 42540090
-p 608 4736286 40655049
-p 607 4736286 40655049
-p 610 4736286 39300412
-p 613 4736286 38409122
-p 614 4736286 37517832
-p 620 6530338 35180385
-p 618 6530338 35180385
-p 619 4736286 35180385
-p 617 4736286 35180385
-p 621 17655148 35180385
-p 626 6530338 33734229
-p 624 6530338 33734229
-p 625 4736286 33734229
-p 623 4736286 33734229
-p 627 6530338 32842939
-p 628 18396802 32842939
-p 630 6530338 31951649
-p 631 12092743 31951649
-p 633 6530338 28386489
-p 634 12092743 28386489
-p 636 6530338 27495199
-p 637 12092743 27495199
-p 639 6530338 25712619
-p 640 6530338 24821329
-p 641 9867781 24821329
-p 646 6530338 22483882
-p 644 6530338 22483882
-p 645 4736286 22483882
-p 643 4736286 22483882
-p 647 14688532 22483882
-p 649 19138456 21592592
-p 651 6530338 20701302
-p 652 15059359 20701302
-p 654 6530338 17136142
-p 655 12092743 17136142
-p 657 6530338 16244852
-p 658 8384473 16244852
-p 660 6530338 14462272
-p 724 4736286 12022365
-p 723 4736286 12022365
-p 726 4736286 10667728
-p 729 4736286 9776438
-p 730 4736286 8885148
-p 772 4736286 6027778
-p 773 4736286 6027778
-l 778 334
-l 779 334
-l 780 334
-l 781 334
-l 782 334
-l 783 334
-l 784 334
-l 785 334
-l 786 334
-l 787 334
-l 788 334
-l 789 334
-l 790 334
-l 791 334
-l 792 334
-l 793 334
-s 9
-p 793 4736286 6027778
-p 790 4736286 50644704
-p 791 4736286 50644704
-p 792 4736286 49530837
-p 782 4736286 50644704
-p 784 4736286 49465056
-p 785 4736286 49465056
-p 786 4736286 48678624
-p 781 4736286 47040224
-p 780 4736286 51336253
-p 779 4736286 51336253
-p 661 4736286 43268483
-p 662 4736286 43268483
-p 668 6620039 36389056
-p 666 6620039 36389056
-p 667 6261229 36389056
-p 665 6261229 36389056
-p 677 6620039 42625056
-p 678 6620039 42625056
-p 669 6620039 36389056
-p 670 6620039 36389056
-p 672 6620039 36389056
-p 673 6620039 36389056
-p 675 6620039 36389056
-p 676 6620039 36389056
-p 674 6620039 36389056
-p 681 6620039 35497766
-p 679 10509803 35497766
-p 682 15409447 36389056
-p 683 15409447 37280346
-p 692 15409447 43268483
-p 693 15409447 43268483
-p 684 15409447 36389056
-p 685 15409447 36389056
-p 687 15409447 36389056
-p 688 15409447 36389056
-p 690 15409447 36389056
-p 691 15409447 36389056
-p 689 15409447 36389056
-p 696 15409447 35497766
-p 694 19256064 35497766
-p 697 24146230 36389056
-p 698 24146230 37280346
-p 707 24146230 43097551
-p 708 24146230 43097551
-p 699 24146230 36389056
-p 700 24146230 36389056
-p 702 24146230 36389056
-p 703 24146230 36389056
-p 705 24146230 36389056
-p 706 24146230 36389056
-p 704 24146230 36389056
-p 711 24146230 35497766
-p 709 27756878 35497766
-p 716 4736286 33763681
-p 717 8585921 33763681
-p 718 8585921 34654971
-p 731 4736286 22947533
-p 732 4736286 22947533
-p 738 7781060 14238183
-p 736 7781060 14238183
-p 737 7422250 14238183
-p 735 7422250 14238183
-p 739 7781060 14238183
-p 741 7781060 14238183
-p 742 7781060 14238183
-p 744 7781060 14238183
-p 745 7781060 14238183
-p 743 7781060 14238183
-p 746 19700503 14238183
-p 748 19700503 14238183
-p 749 19700503 14238183
-p 751 19700503 14238183
-p 752 19700503 14238183
-p 750 19700503 14238183
-p 757 4736286 12691533
-p 754 14051817 12691533
-p 755 17173800 12691533
-p 756 17173800 13582823
-p 788 4736286 6027778
-p 789 4736286 6027778
-l 794 334
-l 795 335
-l 796 336
-l 797 337
-l 798 338
-l 799 339
-l 800 340
-l 801 341
-l 802 344
-l 803 344
-l 804 344
-l 805 345
-l 806 346
-l 807 347
-l 808 348
-l 809 349
-l 810 350
-l 811 353
-l 812 354
-l 813 355
-l 814 356
-l 815 357
-l 816 359
-l 817 359
-l 818 359
-l 819 360
-l 820 361
-l 821 362
-l 822 362
-l 823 362
-l 824 363
-l 825 364
-l 826 365
-l 827 366
-l 828 366
-l 829 366
-l 830 367
-l 831 370
-l 832 372
-l 833 374
-l 834 375
-l 835 376
-l 836 377
-l 837 380
-l 838 380
-l 839 380
-l 840 382
-l 841 382
-l 842 382
-l 843 383
-l 844 383
-l 845 383
-l 846 383
-l 847 383
-l 848 383
-l 849 383
-l 850 383
-l 851 383
-l 852 383
-l 853 383
-l 854 383
-l 855 383
-s 10
-p 855 4736286 6027778
-p 852 4736286 50644704
-p 853 4736286 50644704
-p 854 4736286 49530837
-p 844 4736286 50644704
-p 846 4736286 49465056
-p 847 4736286 49465056
-p 848 4736286 48678624
-p 843 4736286 47040224
-p 764 6530338 46319328
-p 762 6530338 46319328
-p 763 4736286 46319328
-p 761 4736286 46319328
-p 794 6530338 45428038
-p 795 6530338 44536748
-p 796 6530338 43645458
-p 797 6530338 42754168
-p 798 6530338 41862878
-p 799 6530338 40971588
-p 800 6530338 40080298
-p 801 6530338 39189008
-p 805 6530338 37891931
-p 803 6530338 37891931
-p 804 4736286 37891931
-p 802 4736286 37891931
-p 806 6530338 37000641
-p 807 6530338 36109351
-p 808 6530338 35218061
-p 809 6530338 34326771
-p 810 6530338 33435481
-p 811 4736286 32138404
-p 812 4736286 30355824
-p 813 4736286 29464534
-p 814 4736286 28573244
-p 815 4736286 26790664
-p 819 6530338 24602296
-p 817 6530338 24602296
-p 818 4736286 24602296
-p 816 4736286 24602296
-p 820 6530338 23711006
-p 821 6530338 22819716
-p 822 12092743 22819716
-p 824 6530338 21928426
-p 825 6530338 21037136
-p 826 6530338 20145846
-p 827 6530338 19254556
-p 828 9867781 19254556
-p 830 6530338 18363266
-p 831 4736286 17066189
-p 832 4736286 14392319
-p 833 4736286 12609739
-p 834 4736286 10827159
-p 835 4736286 9935869
-p 836 4736286 8153289
-p 850 4736286 6027778
-p 851 4736286 6027778
-l 856 383
-l 857 385
-l 858 388
-l 859 388
-l 860 388
-l 861 389
-l 862 389
-l 863 390
-l 864 391
-l 865 391
-l 866 392
-l 867 393
-l 868 393
-l 869 394
-l 870 395
-l 871 395
-l 872 396
-l 873 399
-l 874 404
-l 875 404
-l 876 404
-l 877 406
-l 878 406
-l 879 406
-l 880 407
-l 881 409
-l 882 411
-l 883 411
-l 884 411
-l 885 411
-l 886 411
-l 887 411
-l 888 411
-l 889 413
-l 890 413
-l 891 413
-l 892 413
-l 893 413
-l 894 413
-l 895 413
-l 896 413
-l 897 413
-l 898 413
-l 899 413
-l 900 413
-l 901 413
-l 902 413
-l 903 414
-l 904 414
-l 905 415
-l 906 415
-l 907 415
-l 908 415
-l 909 415
-l 910 415
-l 911 415
-l 912 415
-l 913 415
-l 914 415
-l 915 415
-l 916 415
-l 917 415
-l 918 416
-l 919 416
-l 920 416
-l 921 416
-l 922 416
-l 923 416
-l 924 416
-l 925 419
-l 926 419
-l 927 419
-l 928 422
-l 929 422
-l 930 422
-l 931 423
-l 932 425
-l 933 428
-l 934 428
-l 935 428
-l 936 429
-l 937 431
-l 938 431
-l 939 431
-l 940 431
-l 941 431
-l 942 431
-l 943 431
-l 944 431
-l 945 431
-l 946 431
-l 947 431
-l 948 431
-l 949 431
-s 11
-p 949 4736286 6027778
-p 946 4736286 50644704
-p 947 4736286 50644704
-p 948 4736286 49530837
-p 938 4736286 50644704
-p 940 4736286 49465056
-p 941 4736286 49465056
-p 942 4736286 48678624
-p 937 4736286 47040224
-p 927 4736286 46319328
-p 925 4736286 46319328
-p 839 4736286 45598432
-p 838 4736286 45598432
-p 841 4736286 43326291
-p 856 4736286 41543711
-p 857 4736286 40652421
-p 861 6530338 38974699
-p 859 6530338 38974699
-p 860 4736286 38974699
-p 858 4736286 38974699
-p 862 6530338 38974699
-p 863 6530338 38083409
-p 864 6530338 37192119
-p 865 6530338 37192119
-p 866 6530338 36300829
-p 867 6530338 35409539
-p 868 6530338 35409539
-p 869 6530338 34518249
-p 870 6530338 33626959
-p 871 6530338 33626959
-p 872 6530338 32735669
-p 873 4736286 31057947
-p 876 4736286 24680075
-p 875 4736286 24680075
-p 878 4736286 23325438
-p 880 4736286 21542858
-p 881 4736286 20651568
-p 930 4736286 12491116
-p 929 4736286 12491116
-p 931 4736286 11136479
-p 932 4736286 10245189
-p 944 4736286 6027778
-p 945 4736286 6027778
-l 950 431
-l 951 433
-l 952 433
-l 953 433
-l 954 435
-l 955 439
-l 956 439
-l 957 439
-l 958 441
-l 959 443
-l 960 443
-l 961 443
-l 962 445
-(.geocubit_manual.bbl
-)
-l 963 447
-l 964 447
-l 965 447
-l 966 447
-l 967 447
-l 968 447
-l 969 447
-l 970 447
-l 971 447
-l 972 447
-l 973 447
-l 974 447
-l 975 447
-l 976 447
-l 977 447
-l 978 447
-s 12
-p 978 4736286 6027778
-p 975 4736286 50644704
-p 976 4736286 50644704
-p 977 4736286 49530837
-p 967 4736286 50644704
-p 969 4736286 49465056
-p 970 4736286 49465056
-p 971 4736286 48678624
-p 966 4736286 47040224
-p 965 4736286 47040224
-p 964 4736286 47040224
-p 882 4736286 47040224
-p 883 4736286 47040224
-p 889 10040492 40096214
-p 887 10040492 40096214
-p 888 9681682 40096214
-p 886 9681682 40096214
-p 898 10040492 46661117
-p 899 10040492 46661117
-p 890 10040492 40096214
-p 891 10040492 40096214
-p 893 10040492 40096214
-p 894 10040492 40096214
-p 896 10040492 40096214
-p 897 10040492 40096214
-p 895 10040492 40096214
-p 902 10040492 39204924
-p 900 13200092 39204924
-p 903 17369571 40096214
-p 904 17369571 40987504
-p 913 17369571 47040224
-p 914 17369571 47040224
-p 905 17369571 40096214
-p 906 17369571 40096214
-p 908 17369571 40096214
-p 909 17369571 40096214
-p 911 17369571 40096214
-p 912 17369571 40096214
-p 910 17369571 40096214
-p 917 17369571 39204924
-p 915 22624651 39204924
-p 922 4736286 37470839
-p 923 8093734 37470839
-p 924 8093734 38362129
-p 935 4736286 32710780
-p 934 4736286 32710780
-p 936 4736286 31356143
-p 950 4736286 30464853
-p 953 4736286 28543431
-p 952 4736286 28543431
-p 954 4736286 27188794
-p 957 4736286 24928003
-p 956 4736286 24928003
-p 958 4736286 23326247
-p 961 4736286 21065456
-p 960 4736286 21065456
-p 962 4736286 17104404
-p 973 4736286 6027778
-p 974 4736286 6027778
-l 979 447
-l 980 447
Modified: seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/doc/geocubit_manual.tex
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/doc/geocubit_manual.tex 2013-09-26 00:49:25 UTC (rev 22849)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/doc/geocubit_manual.tex 2013-09-26 00:57:45 UTC (rev 22850)
@@ -2,82 +2,89 @@
%% LyX 2.0.0beta3 created this file. For more info, see http://www.lyx.org/.
%% Do not edit unless you really know what you are doing.
-\documentclass[a4paper,11pt,english]{article}
-\setlength{\textwidth}{450pt}
-\setlength{\oddsidemargin}{0pt}
-\setlength{\marginparwidth}{0pt}
+\documentclass[a4paper,11pt,oneside,english]{article}
\setcounter{secnumdepth}{3}
\setcounter{tocdepth}{3}
\usepackage[T1]{fontenc}
-\usepackage[latin9]{inputenc}
+\usepackage[utf8]{inputenc}
+\usepackage[margin=2.5cm]{geometry}
\usepackage{babel}
-\usepackage{setspace}
\usepackage{listings}
-\usepackage{graphicx,subfigure}
-\usepackage{pdfsync}
+\usepackage{graphicx,subcaption}
\usepackage{colortbl}
\usepackage[unicode=true,pdfusetitle,
bookmarks=true,bookmarksnumbered=true,bookmarksopen=false,
- breaklinks=false,pdfborder={0 0 0},backref=section,colorlinks=true,pdfpagemode=None]
+ breaklinks=false,pdfborder={0 0 0},backref=section,colorlinks=true]
{hyperref}
-\hypersetup{
- pdftex,pdfstartview=FitH,pdfpagelayout=SinglePage}
+\hypersetup{pdfstartview=FitH,pdfpagelayout=SinglePage}
+\newcommand{\cubit}{\textsc{Cubit}}
+\newcommand{\geocubit}{\textsc{GeoCubit}}
+\newcommand{\specfem}{\textsc{Specfem3D}}
-\newenvironment{lyxcode}
-{\par\begin{list}{}{
-\setlength{\rightmargin}{\leftmargin}
-\setlength{\listparindent}{0pt}% needed for AMS classes
-\setlength{\itemsep}{0pt}
-\setlength{\parsep}{0pt}
-\normalfont\ttfamily}%
-\item[]}
-{\end{list}}
+\definecolor{codebg}{RGB}{238,238,238}
+\definecolor{codeframe}{RGB}{204,204,204}
+\lstset{
+ % The bash style works pretty well for cfg files as well.
+ language=bash,
+ basicstyle=\small\ttfamily,
+ breakatwhitespace=true,
+ breaklines=true,
+ escapechar=\@,
+ backgroundcolor=\color{codebg},
+ frame=single,
+ framesep=10pt,
+ rulecolor=\color{codeframe}
+}
\begin{document}
-\title{GEOCUBIT USER MANUAL}
+\title{\geocubit\ User Manual}
\author{Emanuele Casarotti}
\maketitle
\tableofcontents
-\newpage{}
-\section{WHY}
+\clearpage
+\section{Why?}
-GEOCUBIT is a python library wrapping around the \href{http://cubit.sandia.gov} {CUBIT} Python Interface and it aims to facilitate the meshing
-process in some common problems in seismic wave propagation.\\
-In particular, it is focused on the meshing requests of SPECFEM3D and it is helpful for some tedious tasks as:\\
-
+\geocubit\ is a Python library wrapping around the \href{http://cubit.sandia.gov} {\cubit}
+Python Interface and it aims to facilitate the meshing
+process in some common problems in seismic wave propagation.
+In particular, it is focused on the meshing requests of \specfem\ and it is helpful for such tedious tasks as:
+%
\begin{itemize}
\item Creation of geophysical surfaces and volumes (ex. topography).
\item Mesh of layered volumes with hexahedral.
\item Creation of an anisotropic mesh suitable for cases where some alluvial
basin (or slow velocity zones) are present.
-\item It can be used as serial or parallel process. The parallel meshing
+\item It can be used as a serial or parallel process. The parallel meshing
capabilities are fundamental for large geophysical problems (ex. mesh
of Southern California using SRTM topography).
\end{itemize}
-GEOCUBIT can be use inside the graphical interface of \href{http://cubit.sandia.gov} {CUBIT} (i.e.
-as python object in the script tab) or as unix command. We refer to the \href{http://cubit.sandia.gov} {CUBIT} help appendix for more information about the python interface.
+\geocubit\ can be used inside the graphical interface of \cubit\ (i.e.
+as a Python object in the script tab) or as unix command. We refer to the
+\cubit\ help appendix for more information about the Python interface.
-\newpage{}
-\section{FIRST STEPS}
+\clearpage
+\section{First Steps}
\subsection{Requirements}
-The minimum requirements for using GEOCUBIT are:
-
+The minimum requirements for using \geocubit\ are:
+%
\begin{itemize}
-\item \href{http://cubit.sandia.gov} {CUBIT 12.2}
-\item \href{http://downloads.sourceforge.net/numpy} {numpy 1.0+}
-\item \href{http://www.python.org} {python 2.5}
- (strictly! it depends on the cubit library that refers to this version of python)
+\item \href{http://cubit.sandia.gov} {\cubit\ 12.2}
+\item \href{http://downloads.sourceforge.net/numpy} {NumPy 1.0+}
+\item \href{http://www.python.org} {Python 2.5}
+ (strictly! it depends on the \cubit\ library that refers to this version of Python)
\end{itemize}
-\raggedright{and for using the parallel meshing capabilities:}
+%
+and for using the parallel meshing capabilities:
+%
\begin{itemize}
\item \href{http://downloads.sourceforge.net/pympi/pyMPI-2.4b4.tar.gz} {pyMPI}
\end{itemize}
@@ -85,386 +92,398 @@
\subsection{Installation}
For installing, download the code and type in the GEOCUBIT directory:
-\begin{lyxcode}
-python2.5 setup.py install~
-\end{lyxcode}
+%
+\begin{lstlisting}
+python2.5 setup.py install
+\end{lstlisting}
+%
+and check that the following variables are set:
+%
+\begin{lstlisting}
+CUBITDIR="/usr/local/CUBIT"
+CUBITLIB="$CUBITDIR/bin:$CUBITDIR/structure:$CUBITDIR/components"
+PYTHONPATH="$CUBITDIR/components/cubit:$CUBITDIR/structure:$CUBITDIR/bin"
+LD_LIBRARY_PATH="$CUBITDIR/bin"
+PATH="$CUBITDIR/bin:$CUBITDIR/components/cubit"
+\end{lstlisting}
-Check that the following variables are set:
-\begin{lyxcode}
-CUBITDIR=/usr/local/CUBIT\\
-CUBITLIB=\$CUBITDIR/bin:\$CUBITDIR/structure:\$CUBITDIR/components\\
-PYTHONPATH=\$CUBITDIR/components/cubit:\$CUBITDIR/structure:\$CUBITDIR/bin\\
-LD\_LIBRARY\_PATH=\$CUBITDIR/bin\\
-PATH=\$CUBITDIR/bin:\$CUBITDIR/components/cubit \\
-\end{lyxcode}
+\clearpage
+\section{Using \geocubit\ from the Command Line}
+\subsection{Utilities}
-\newpage{}
-\section{USING GEOCUBIT AT COMMAND LINE:}
-
-\subsection{UTILITIES}
\begin{itemize}
-\item checking the configuration of the libraries and dependencies:\\
-\begin{lyxcode}
+\item Checking the configuration of the libraries and dependencies:
+\begin{lstlisting}
GEOCUBIT.py --chklib
-\end{lyxcode}
+\end{lstlisting}
-\item checking the parameter file:\\
-\begin{lyxcode}
+\item Checking the parameter file:
+\begin{lstlisting}
GEOCUBIT.py --chkcfg --cfg=[file]
-\end{lyxcode}
+\end{lstlisting}
\end{itemize}
-\subsection{CREATE GEOMETRY}
-SURFACES
+\subsection{Create Geometry}
+\subsubsection{Surfaces}
+
\begin{itemize}
-\item creating acis surfaces using a parameter file:\\
-\begin{lyxcode}
- GEOCUBIT.py --build\_surface --cfg=[filename]
-\end{lyxcode}
+\item Creating acis surfaces using a parameter file:
+\begin{lstlisting}
+ GEOCUBIT.py --build_surface --cfg=[filename]
+\end{lstlisting}
\end{itemize}
\begin{itemize}
-\item creating a surface from regular ascii grid or ascii lines that define a "skin":\\
-\begin{lyxcode}
+\item Creating a surface from regular ascii grid or ascii lines that define a "skin":
+\begin{lstlisting}
GEOCUBIT.py --surface=[file] (--regulargrid=[opts] | --skin=[opts])
-\end{lyxcode}
+\end{lstlisting}
\end{itemize}
-VOLUMES
+\subsubsection{Volumes}
+
\begin{itemize}
-\item **serial** command: creating cubit volumes using a parameter file:\\
-\begin{lyxcode}
+\item \emph{Serial} command: creating \cubit\ volumes using a parameter file:
+\begin{lstlisting}
GEOCUBIT.py --surface=[file] (--regulargrid=[opts] | --skin=[opts])
-\end{lyxcode}
+\end{lstlisting}
\end{itemize}
\begin{itemize}
-\item **parallel** command: creating cubit volumes using a parameter file:\\
-\begin{lyxcode}
-mpirun -n [numproc] pyMPI GEOCUBIT.py --build\_volume --cfg=[file]
-\end{lyxcode}
+\item \emph{Parallel} command: creating \cubit\ volumes using a parameter file:
+\begin{lstlisting}
+mpirun -n [numproc] pyMPI GEOCUBIT.py --build_volume --cfg=[file]
+\end{lstlisting}
-\begin{figure}\begin{center}
- \subfigure[]{\includegraphics[scale=.5]{media/geocubitmanual_grid2.jpg}
- \label{fig:vol}}%
- \subfigure[]{\includegraphics[scale=.9]{media/geocubitmanual_grid.jpg}
- \label{fig:table}}%
-\caption{a) Volume divided in slices. b) Parallel map of the IDs for volume in Figure~\ref{fig:vol} divided in (16 slices, 4 along x-axis and 4 along y-axis).}
-\label{fig:grid}
-\end{center}
+\begin{figure}[htbp]
+ \centering
+ \begin{subfigure}{0.45\textwidth}
+ \includegraphics[width=\textwidth,keepaspectratio]{media/geocubitmanual_grid2.jpg}
+ \caption{Volume divided in slices.}
+ \label{fig:vol}
+ \end{subfigure}
+ \begin{subfigure}{0.45\textwidth}
+ \includegraphics[width=\textwidth,keepaspectratio]{media/geocubitmanual_grid.jpg}
+ \caption{Parallel map of the IDs for volume in \autoref{fig:vol} divided in 16 slices (4 along x-axis and 4 along y-axis).}
+ \label{fig:table}
+ \end{subfigure}
+ \caption{Volume slicing.}
+ \label{fig:grid}
\end{figure}
-In this parallel application, the volume is separated in N slices, each one is assigned at a single process and one file for each slice is created (see Figure~\ref{fig:grid} as example).
+In this parallel application, the volume is separated in $N$ slices, each one is assigned at a single process and one file for each slice is created (see \autoref{fig:grid} as example).
\end{itemize}
-\subsection{MESHING}
-\begin{itemize}
-\item **serial** command: building a volume and meshing it.\\
-\begin{lyxcode}
-GEOCUBIT.py --build\_volume --mesh --cfg=[file] (--id\_proc=n)
-\end{lyxcode}
-\end{itemize}
+\subsection{Meshing}
\begin{itemize}
-\item **serial** command : meshing a volume\\
-\begin{lyxcode}
-GEOCUBIT.py --mesh --cfg=[file] (--id\_proc=[n, def=0])
-\end{lyxcode}
-\footnotesize note: without the --build\_volume flag the script recall an old "geometry\_vol\_[id\_proc].cub" file\\
-\end{itemize}
+\item \emph{Serial} command: building a volume and meshing it.
+\begin{lstlisting}
+GEOCUBIT.py --build_volume --mesh --cfg=[file] (--id_proc=n)
+\end{lstlisting}
-\begin{itemize}
-\item **parallel** command: meshing a volume from a parameter file\\
-\begin{lyxcode}
-mpirun -n [nproc] pyMPI GEOCUBIT.py --build\_volume --mesh --cfg=[file]
-\end{lyxcode}
+\item \emph{Serial} command: meshing a volume
+\begin{lstlisting}
+GEOCUBIT.py --mesh --cfg=[file] (--id_proc=[n, def=0])
+\end{lstlisting}
+
+{\footnotesize note: without the \lstinline|--build_volume| flag the script recall an old "\lstinline|geometry_vol_[id_proc].cub|" file}
+
+\item \emph{Parallel} command: meshing a volume from a parameter file
+\begin{lstlisting}
+mpirun -n [nproc] pyMPI GEOCUBIT.py --build_volume --mesh --cfg=[file]
+\end{lstlisting}
\end{itemize}
-\subsection{FINALIZING AND EXPORTING}
+\subsection{Finalizing and Exporting}
-It is possible to collect, merge, set the absorbing boundary conditions and to export the resulting mesh in a format readable by SPECFEM3D.
+It is possible to collect, merge, set the absorbing boundary conditions and to export the resulting mesh in a format readable by \specfem.
\begin{itemize}
\item
-Collecting some cubit files, setting the absorbing boundary conditions and merging in a single free mesh cubitfile
-\begin{lyxcode}
+Collecting some \cubit\ files, setting the absorbing boundary conditions and merging in a single free mesh cubitfile
+\begin{lstlisting}
GEOCUBIT.py --collect --merge --meshfiles=[files] --cpux=N --cpuy=N (--rangecpux=[cpuxmin,cpuxmax] --rangecpuy=[cpuymin,cpuymax])
-\end{lyxcode}
-\texttt{cpuy, cpux} set the number of cpus used for creating the mesh files; \texttt{rangecpux} and \texttt{rangecpuy} set the range of slices that are used in the mesh.\\
-Following the example in Figure~\ref{fig:grid}, the parameters \texttt{--cpux=4 --cpuy=4 --rangecpux=1,2 --rangecpuy=2,3} select the slices with id 10 11 14 15. Only these slices are going to be collect and merged with the appropriated absorbing boundary conditions.
+\end{lstlisting}
-\item Collecting a single free mesh cubitfile and refine the hex inside some curves (ex. alluvial basin, Figure~\ref{fig:refcurve})
-\begin{lyxcode}
-GEOCUBIT.py --collect --meshfiles=[list of files] --curverefining=[file]\\
-\end{lyxcode}
-\footnotesize note: the curves must be stored in acis format (sat) and must be closed.
+\lstinline{cpux, cpuy} set the number of cpus used for creating the mesh files; \lstinline{rangecpux} and \lstinline{rangecpuy} set the range of slices that are used in the mesh.
-\end{itemize}
+Following the example in \autoref{fig:grid}, the parameters \lstinline{--cpux=4 --cpuy=4 --rangecpux=1,2 --rangecpuy=2,3} select the slices with id $10, 11, 14, 15$. Only these slices are going to be collected and merged with the appropriate absorbing boundary conditions.
-\begin{figure}\center
-\includegraphics[scale=.37]{media/refcurve_1.png}%
-\includegraphics[scale=.303]{media/refcurve_2.png}%
-\caption{}
+\item Collecting a single free mesh cubitfile and refine the hex inside some curves (ex. alluvial basin, \autoref{fig:refcurve})
+\begin{lstlisting}
+GEOCUBIT.py --collect --meshfiles=[list of files] --curverefining=[file]
+\end{lstlisting}
+{\footnotesize note: the curves must be stored in acis format (sat) and must be closed.}
+
+\begin{figure}[htbp]
+\center
+\includegraphics[height=0.4\textheight,keepaspectratio]{media/refcurve_1.png}%
+\includegraphics[height=0.4\textheight,keepaspectratio]{media/refcurve_2.png}%
+\caption{An alluvial basin example.}
+
\label{fig:refcurve}
\end{figure}
-\begin{itemize}
-\item Exporting a cubit mesh file in a SPECFEM3D mesh
-\begin{lyxcode}
-GEOCUBIT.py --export2SPECFEM3D=[output directory] --meshfiles=[filename] (--listblock=[list of cubit blocks] --listflag=[list of specfem3d flags])
-\end{lyxcode}
-\footnotesize If required, it is possible to assign personalized flags to each block using --listblock and the correspondent listflag (for example: --listblock=3,4 --listflag=1,-1)
+\item Exporting a \cubit\ mesh file to a \specfem\ mesh
+\begin{lstlisting}
+GEOCUBIT.py --export2SPECFEM3D=[output directory] --meshfiles=[filename] (--listblock=[list of @\cubit@ blocks] --listflag=[list of specfem3d flags])
+\end{lstlisting}
+
+If required, it is possible to assign personalized flags to each block using \lstinline{--listblock} and the corresponding listflag (for example: \lstinline{--listblock=3,4 --listflag=1,-1})
\end{itemize}
-\newpage{}
-\section{GEOCUBIT and CUBIT GRAPHICAL INTERFACE}
-In the python script tab of CUBIT GUI, type:
+\clearpage
+\section{\geocubit\ and \cubit\ Graphical Interface}
-\begin{lyxcode}
->f="volume.cfg" $\to$ \textrm{\small{\textit{store the name of the configuration file, see next section}}}\\
-\end{lyxcode}
+In the Python script tab of \cubit\ GUI, type:
-\begin{lyxcode}
->from geocubitlib import volumes $\to$ \textrm{\small{\textit{load the geocubit modules}}}\\
->from geocubitlib import mesh\_volume\\
->from geocubitlib import exportlib\\
-\end{lyxcode}
+\begin{lstlisting}
+>f="volume.cfg" # store the name of the configuration file, see next section
-\begin{lyxcode}
-volumes.volumes(f) $\to$ \textrm{\small{\textit{create the volumes, the parameters are stored in the cfg file}}}\\
-mesh\_volume.mesh(f) $\to$ \textrm{\small{\textit{mesh the volume}}}\\
-\end{lyxcode}\begin{lyxcode}
-exportlib.collect() $\to$ \textrm{\small{\textit{set the boundary conditions, see SPECFEM3D manual}}}\\
-exportlib.e2SEM(outdir="./output") $\to$ \textrm{\small{\textit{save the mesh in the SPECFEM3D format}}}\\
-\end{lyxcode}
-See the media/iterative.mov for a live example.
+>from geocubitlib import volumes # load the @\geocubit@ modules
+>from geocubitlib import mesh_volume
+>from geocubitlib import exportlib
+volumes.volumes(f) # create the volumes, the parameters are stored in the cfg file
+mesh_volume.mesh(f) # mesh the volume
-\section{EXAMPLES of CONFIGURATION FILES}
-In the current section a few configuration files will be described. Please, see the files in the examples directory of the GEOCUBIT distribution for reference.\\
+exportlib.collect() # set the boundary conditions, see @\specfem@ manual
+exportlib.e2SEM(outdir="./output") # save the mesh in the @\specfem@ format
+\end{lstlisting}
-\subsection{GENERAL OPTIONS}
+See the \lstinline|media/iterative.mov| for a live example.
+\section{Examples of Configuration Files}
+
+In the current section a few configuration files will be described. Please see the files in the examples directory of the \geocubit\ distribution for reference.
+
+\subsection{General Options}
+
A cfg file has a format similar to ini windows files.
-They are divided in \texttt{[sections]}. There is no order in the position or in the name of the sections.
+They are divided in \lstinline{[sections]}. There is no order in the position or in the name of the sections.
-Generally, the parameters that control the general options are in the first section, called \texttt{[cubit.options]} and \texttt{[simulation.cpu\_parameters]}. For example:
+Generally, the parameters that control the general options are in the first section, called \lstinline{[cubit.options]} and \lstinline{[simulation.cpu_parameters]}. For example:
-\begin{lyxcode}
-[cubit.options]\\
-cubit\_info=off $\to$ \textrm{\small{\textit{turn on/off the information of the cubit command}}}\\
-echo\_info=off $\to$ \textrm{\small{\textit{turn on/off the echo of the cubit command}}}\\
-cubit\_info=off $\to$ \textrm{\small{\textit{turn on/off the cubit journaling}}}\\
-echo\_info=off $\to$ \textrm{\small{\textit{turn on/off the cubit error journaling}}}\\
-working\_dir=tmp $\to$ \textrm{\small{\textit{set the working directory}}}\\
-output\_dir=output $\to$ \textrm{\small{\textit{set the output directory}}}\\
-save\_geometry\_cubit = True $\to$ \textrm{\small{\textit{true if it saves the geometry files}}}\\
-save\_surface\_cubit = False $\to$ \textrm{\small{\textit{true if it saves the surfaces in a cubit files}}}\\
-export\_exodus\_mesh = True $\to$ \textrm{\small{\textit{true if it saves the mesh also in exodus format. The default is the cubit format (cub), that contains both the geometry and the meshing information. In case of complex geometry, the cub file could be enormous and a more light exodus file become important.}}}\\
-\end{lyxcode}
-\begin{lyxcode}
-[simulation.cpu\_parameters]\\
-nodes = 1 $\to$ \textrm{\small{\textit{number of nodes/process}}}\\
+\begin{lstlisting}
+[cubit.options]
+cubit_info=off # turn on/off the information of the @\cubit@ command
+echo_info=off # turn on/off the echo of the @\cubit@ command
+cubit_info=off # turn on/off the @\cubit@ journaling
+echo_info=off # turn on/off the @\cubit@ error journaling
+working_dir=tmp # set the working directory
+output_dir=output # set the output directory
+save_geometry_cubit = True # true if it saves the geometry files
+save_surface_cubit = False # true if it saves the surfaces in a @\cubit@ files
+export_exodus_mesh = True # true if it saves the mesh also in exodus format. The default is the @\cubit@ format (cub), that contains both the geometry and the meshing information. In case of complex geometry, the cub file could be enormous and a more light exodus file become important.
-\end{lyxcode}
+[simulation.cpu_parameters]
+nodes = 1 # number of nodes/process
+\end{lstlisting}
-\footnotesize note: Usually the *\_info options are turn off by default for improving the performances.\normalsize\\
+{\footnotesize note: Usually the \lstinline|*_info| options are turned off by default for improving the performance.}
-\subsection{CREATING A TOPOGRAPHIC SURFACE}
-\footnotesize See stromboli.cfg, execute with: \verb!GEOCUBIT.py --build\_surface --cfg=./example/stromboli.cfg. !\normalsize\\
-~\\
-GEOCUBIT is able to create a topographic surface based upon the CUBIT command \verb!create skin curve! and \verb!create surface net! (Figure~\ref{fig:surfacetype}). See the CUBIT manual for details.
+\subsection{Creating a Topographic Surface}
-\begin{lyxcode}
+{\footnotesize See stromboli.cfg, execute with: \lstinline|GEOCUBIT.py --build_surface --cfg=./example/stromboli.cfg|.}
+
+\geocubit\ is able to create a topographic surface based upon the \cubit\ command \lstinline|create skin curve| and \lstinline|create surface net| (\autoref{fig:surfacetype}). See the \cubit\ manual for details.
+
+\begin{lstlisting}
[geometry.surfaces]
-nsurf = 2 $\to$ \textrm{\small{\textit{number of surfaces that will be created}}}\\
-\end{lyxcode}
-\begin{lyxcode}
-[surface\textbf{1}.parameters]\\
-name=example/data/stromboli.xyz $\to$ \textrm{\small{\textit{the name of the file defining the surface}}}\\
-surf\_type=skin $\to$ \textrm{\small{\textit{the type of surface:} \texttt{skin} \textit{, the file is a sequence of parallel points that span the surface in the order describer by the following direction parameters} \texttt{regular\_grid} \emph{the points are structured distritibuted. see Cubit Manual for details(skin surface and net surface) }}}\\
-unit\_surf= utm $\to$ \textrm{\small{\textit{unit of the points,} \texttt{utm} \emph{or} \texttt{geo} \emph{(geographical) coordinates}}}\\
-directionx = 0 $\to$ \textrm{\small{\textit{if} \texttt{surf\_type=skin}\emph{:} \texttt{directionx} \emph{and} \texttt{directiony} \emph{define how the points span the surface}}}\\
-directiony = 1\\
-step = 1 $\to$ \textrm{\small{\textit{the script creates a vertex each} \texttt{step} \emph{points: in this case} \texttt{step=1} \textit{means a vertex for all the points}}}\\
-\end{lyxcode}
+nsurf = 2 # number of surfaces that will be created
-\begin{lyxcode}
-[surface2.parameters] $\to$ \textrm{\small{\emph{note the different number of surface in the section name}}}
-name=./examples/surfaces/topo.dat $\to$ \textrm{\small{\textit{the name of the file defining the surface}}}\\
-surf\_type=regular\_grid $\to$ \textrm{\small{\textit{the type of surface:} \texttt{skin} \textit{, the file is a sequence of parallel points that span the surface in the order describer by the following direction parameters} \texttt{regular\_grid} \emph{the points are structured distritibuted. See see Cubit Manual for details (skin surface and net surface) }}}\\
-unit\_surf= geo $\to$ \textrm{\small{\textit{unit of the points,} \texttt{utm} \emph{or} \texttt{geo} \emph{(geographical) coordinates}}}\\
-nx=5 $\to$ \textrm{\small{\textit{In case of} \texttt{regular\_grid}: \texttt{nx} \emph{and} \texttt{ny} \emph{are the number of points along x and y direction.}}}\\
+[surface1.parameters]
+name=example/data/stromboli.xyz # the name of the file defining the surface
+surf_type=skin # the type of surface: @\texttt{skin}@: the file is a sequence of parallel points that span the surface in the order described by the following direction parameters; @\texttt{regular\_grid}@: the points are structured distributed. See @\cubit@ Manual for details (skin surface and net surface)
+unit_surf = utm # unit of the points, @\texttt{utm}@ or @\texttt{geo}@ (geographical) coordinates
+directionx = 0 # if surf_type=skin: @\texttt{directionx}@ and @\texttt{directiony}@ define how the points span the surface
+directiony = 1
+step = 1 # the script creates a vertex each @\texttt{step}@ points: in this case @\texttt{step=1}@ means a vertex for all the points
+
+[surface2.parameters] # note the different number of surface in the section name
+name=./examples/surfaces/topo.dat # the name of the file defining the surface
+surf_type=regular_grid # the type of surface: @\texttt{skin}@: the file is a sequence of parallel points that span the surface in the order described by the following direction parameters; @\texttt{regular\_grid}@: the points are structured distributed. See see @\cubit@ Manual for details (skin surface and net surface)
+unit_surf= geo # unit of the points, @\texttt{utm}@ or @\texttt{geo}@ (geographical) coordinates
+nx=5 # In case of @\texttt{regular\_grid}@: @\texttt{nx}@ and @\texttt{ny}@ are the number of points along x and y direction.
ny=5
-\end{lyxcode}
+\end{lstlisting}
-\begin{figure}\begin{center}
- \subfigure[]{\includegraphics[scale=.4]{media/skinsurface.png}
- \label{fig:skin}}%
- \subfigure[]{\includegraphics[scale=.4]{media/netsurface.png}
- \label{fig:net}}%
- \subfigure[]{\includegraphics[scale=.25]{media/strombolisurface.png}
- \label{fig:stromboli}}%
-\caption{GEOCUBIT creates a topographic surface in 2 ways: a) skin topography, \texttt{directionx=0} and \texttt{directiony=1}, b) topography from a net of structured distritibuted points, \texttt{nx=4} and \texttt{ny=4}, c) Topography of Stromboli volcano (Italy) using the net surface method.}
-\label{fig:surfacetype}
-\end{center}
+\begin{figure}[htbp]
+ \centering
+ \begin{subfigure}{0.4\textwidth}
+ \includegraphics[width=\textwidth,keepaspectratio]{media/skinsurface.png}
+ \caption{Skin topography, \texttt{directionx=0} and \texttt{directiony=1}}
+ \label{fig:skin}
+ \end{subfigure}
+ \begin{subfigure}{0.4\textwidth}
+ \includegraphics[width=\textwidth,keepaspectratio]{media/netsurface.png}
+ \caption{Topography from a net of structured distributed points, \texttt{nx=4} and \texttt{ny=4}}
+ \label{fig:net}
+ \end{subfigure}
+
+ \begin{subfigure}{0.45\textwidth}
+ \includegraphics[width=\textwidth,keepaspectratio]{media/strombolisurface.png}
+ \caption{Topography of Stromboli volcano (Italy) using the net surface method.}
+ \label{fig:stromboli}
+ \end{subfigure}
+ \caption{Two methods of creating a topographic surface in \geocubit.}
+ \label{fig:surfacetype}
\end{figure}
-\subsection{CREATING and MESHING A VOLUME}
-\footnotesize see volume.cfg, for executioning: \verb!GEOCUBIT.py --build\_volume --mesh --cfg=./example/volume.cfg !\normalsize\\
-~\\
-The example is for a volume with topography with dimension of the hexahedra increasing with depth by means of a refinement layer (Figure~\ref{fig:volume}).
+\subsection{Creating and Meshing a Volume}
-\begin{figure}\begin{center}
-\includegraphics[scale=.4]{media/geometryvolume.png}%
-\includegraphics[scale=.377]{media/meshvolume.png}%
+{\footnotesize see volume.cfg, for executioning: \lstinline|GEOCUBIT.py --build_volume --mesh --cfg=./example/volume.cfg|}
+
+The example is for a volume with topography with dimension of the hexahedra increasing with depth by means of a refinement layer (\autoref{fig:volume}).
+
+\begin{figure}[htbp]
+\centering
+\includegraphics[height=0.24\textheight,keepaspectratio]{media/geometryvolume.png}%
+\includegraphics[height=0.24\textheight,keepaspectratio]{media/meshvolume.png}%
\caption{Mesh of a simple volume}
\label{fig:volume}
-\end{center}
\end{figure}
+\begin{lstlisting}
+[geometry.volumes]
+volume_type = layercake_volume_ascii_regulargrid_regularmap
+latitude_min = 13.879
+latitude_max = 14.279
+longitude_min = 40.619
+longitude_max = 40.969
+nx = 5
+ny = 5
+unit = geo
+[geometry.volumes.layercake]
+nz = 2
+bottomflat = True
+depth_bottom = -7000
+filename = ./example/topo.dat
+geometry_format=ascii
+\end{lstlisting}
-\begin{lyxcode}
-[geometry.volumes]\\
-volume\_type = layercake\_volume\_ascii\_regulargrid\_regularmap\\
-latitude\_min = 13.879\\
-latitude\_max = 14.279\\
-longitude\_min = 40.619\\
-longitude\_max = 40.969\\
-nx = 5 \\
-ny = 5\\
-unit = geo\\
-\end{lyxcode}
+\begin{description}
+\item[nz] is the number of the horizontal surfaces in the volume (in this case: topography and bottom).
+\item[bottomflat] is \lstinline|True| if the bottom is flat.
+\item[depth\_bottom] is the depth of the bottom.
+\item[filename] is a list of files. Each one defines a surface (in this case there is only the topography file in a \lstinline|regular_grid| format).
+\item[geometry\_format] is set to \lstinline|ascii| since the definition of the surfaces comes from ASCII files (structured xyz points).
+\end{description}
-\begin{lyxcode}
-[geometry.volumes.layercake]\\
-nz = 2\\
-bottomflat = True\\
-depth\_bottom = -7000\\
-filename = ./example/topo.dat\\
-geometry\_format=ascii\\
-\end{lyxcode}
+\begin{lstlisting}
+[meshing]
+map_meshing_type=regularmap
+iv_interval=5, # if only one value is present, append the comma
+size=2000
+or_mesh_scheme=map
+ntripl=1
+tripl=2, # if only one value is present, append the comma
+smoothing=False
+\end{lstlisting}
-\texttt{nz} is the number of the horizontal surfaces in the volume (in this case: topography and bottom).\\
-\texttt{bottomflat} is True if the bottom is flat.\\
-\texttt{depth\_bottom} is the depth of the bottom.\\
-\texttt{filename} is a list of files. Each one defines a surface (in this case there is only the topography file in a regular\_grid format).\\
-\texttt{geometry\_format} is set to ascii since the definition of the surfaces comes from ascii files (structured xyz points). \\
+The \lstinline|[meshing]| section contains the parameters requested for the meshing process.
-\begin{lyxcode}
-[meshing]\\
-map\_meshing\_type=regularmap\\
-iv\_interval=5, $\to$ \textrm{\small{\textit{if only one value is present, append the comma}}}\\
-size=2000\\
-or\_mesh\_scheme=map\\
-ntripl=1\\
-tripl=2, $\to$ \textrm{\small{\textit{if only one value is present, append the comma}}}\\
-smoothing=False\\
-\end{lyxcode}
+\begin{description}
+\item[map\_meshing\_type] sets the meshing scheme and is \lstinline|regularmap| by default (other schemes are in preliminary phase).
+\item[iv\_interval] sets the number of "horizontal" hex sheets for each layer.
+\item[size] is the dimension of hexahedral (horizontally).
+\item[or\_mesh\_scheme] is the meshing scheme for the topography (\lstinline|map| or \lstinline|pave| are possible, see the \cubit\ manual for more information).
+\item[ntripl] is the number of tripling layer in the mesh (in this case 1).
+\item[tripl] means in this case that the refinement layer is located at the second surface (the topography). The surfaces are ordered from the bottom (surface 1) to the top (surface 2).
+\item[smoothing] performs the smoothing command in \cubit.
+\end{description}
-The \texttt{meshing} section contains the parameters request for the meshing process.
-\texttt{map\_meshing\_type} set the meshing scheme and is regularmap by default (other schemes are in preliminary phase).\\
-\texttt{iv\_interval} set the number of "horizontal" hex sheets for each layer.
-\texttt{size} is the dimension of hexahedral (horizontally)\\
-\texttt{or\_mesh\_scheme} is the meshing scheme for the topography (\texttt{map} or \texttt{pave} are possible, see the CUBIT manual for more information).\\
-\texttt{ntripl} is the number of tripling layer in the mesh (in this case 1).\\
-\texttt{tripl} means in this case that the the refinement layer is located at the second surface (the topography). The surfaces are ordered from the bottom (surface 1) to the top (surface 2).\\
-\texttt{smoothing} performes the smoothing command in Cubit.\\
+\subsection{Layered Volume and Mesh for Central Italy (in Parallel)}
-\subsection{LAYERED VOLUME and MESH FOR CENTRAL ITALY (PARALLEL)}
+In the example abruzzo.cfg, the layers are 2. For execution, type: \lstinline|mpirun -n 150 GEOCUBIT.py --build_volume --mesh --cfg=./example/abruzzo.cfg|.
-\footnotesize In the example abruzzo.cfg, the layers are 2. For execution type: \texttt{mpirun -n 150 GEOCUBIT.py --build\_volume --mesh --cfg=./example/abruzzo.cfg}.\normalsize\\
-~\\
+Comparing the previous cfg files, there are few modifications:
+%
+\begin{lstlisting}
+. . .
+nz = 3
+. . .
+filename = example/data/moho_int.xyz,example/data/topo.xyz
+. . .
+iv_interval=8,8 (one value for each layer)
+. . .
+refinement_depth=1,
+\end{lstlisting}
-Comparing the previous cfg files, there are few modification:
-
-
-\begin{lyxcode}
-$\cdots $\\
-nz = 3 \\
-$\cdots $\\
-filename = example/data/moho\_int.xyz,example/data/topo.xyz\\
-$\cdots $\\
-iv\_interval=8,8 (one value for each layer)\\
-$\cdots $\\
-refinement\_depth=1,\\
-\end{lyxcode}
-
-The volume has \texttt{nz=3} interfaces: bottom, moho, and topography. The bottom is flat and his z-coordinate position is set by \texttt{depth\_bottom}.
+The volume has \lstinline|nz=3| interfaces: bottom, moho, and topography. The bottom is flat and its z-coordinate position is set by \lstinline|depth_bottom|.
The name of the files that storage the data for the other interfaces are listed in \texttt{filename}.
-The interfaces defines 2 geological layers. Each layer has \texttt{iv\_interval=8} "horizontal" hex sheets. There is only a refinement and it is set in the \texttt{refinement\_depth=1} hex sheet, i.e just below the topography (\texttt{refinement\_depth=8} means just below the moho).
+The interfaces defines 2 geological layers. Each layer has \lstinline|iv_interval=8| "horizontal" hex sheets. There is only a refinement and it is set in the \lstinline|refinement_depth=1| hex sheet, i.e just below the topography (\lstinline|refinement_depth=8| means just below the moho).
-\subsection{CREATION OF A REGULAR MESH}
+\subsection{Creation of a Regular Mesh}
-\footnotesize In the example grid.cfg, the layers are 3. For execution type: \texttt{GEOCUBIT.py --build\_volume --mesh --cfg=./example/grid.cfg}.\normalsize\\
-~\\
+In the example grid.cfg, the layers are 3. For execution, type: \lstinline|GEOCUBIT.py --build_volume --mesh --cfg=./example/grid.cfg|.
-In this example we create a simple layercake box \texttt{geometry\_format=regmesh} with 3 layers at depth defined by \texttt{zdepth=-7000,-3000,-600,0}. The initial mesh has hex with horizontal size \texttt{size=2000} and the numbers of vertical hex sheets is iv\_interval=3,1,1 respectively for the bottom, middle and top layer. We include a refinement layer (\texttt{ntripl=1}) at \texttt{tripl=2,} interface (the second from the bottom). Since the refinement occurs on the vertical direction, the vertical dimensions of the hexes in the top layer is too small and the quality of the mesh decreases. \texttt{coarsening\_top\_layer=True} remesh the top layer and the number of vertical hex sheet in the vertical is defined by \texttt{actual\_vertical\_interval\_top\_layer=1}. (see Figure~\ref{fig:regrid})
+In this example we create a simple layercake box \lstinline|geometry_format=regmesh| with 3 layers at depth defined by \lstinline|zdepth=-7000,-3000,-600,0|. The initial mesh has hex with horizontal size \lstinline|size=2000| and the numbers of vertical hex sheets is \lstinline|iv_interval=3,1,1| respectively for the bottom, middle and top layer. We include a refinement layer (\lstinline|ntripl=1|) at \lstinline|tripl=2,| interface (the second from the bottom). Since the refinement occurs on the vertical direction, the vertical dimensions of the hexes in the top layer is too small and the quality of the mesh decreases. \lstinline|coarsening_top_layer=True| remesh the top layer and the number of vertical hex sheet in the vertical is defined by \lstinline|actual_vertical_interval_top_layer=1|. (see \autoref{fig:regrid})
-\begin{figure}\begin{center}
- \subfigure[]{\includegraphics[scale=.2]{media/regrid.png}
- \label{fig:regrid_nocoarse}}%
- \subfigure[]{\includegraphics[scale=.47]{media/regrid_coarse.png}
- \label{fig:regrid_coarse}}%
-\caption{Regular mesh with one tripling layer, a) mesh with coarsening\_top\_layer=False, b) mesh with \texttt{coarsening\_top\_layer=True}, the top layer is remeshed with \texttt{actual\_vertical\_interval\_top\_layer=1} vertical hex sheet.}
-\label{fig:regrid}
-\end{center}
+\begin{figure}[htbp]
+ \centering
+ \begin{subfigure}{0.45\textwidth}
+ \includegraphics[width=\textwidth,keepaspectratio]{media/regrid.png}
+ \caption{Mesh with \lstinline|coarsening_top_layer=False|}
+ \label{fig:regrid_nocoarse}
+ \end{subfigure}
+ \begin{subfigure}{0.45\textwidth}
+ \includegraphics[width=\textwidth,keepaspectratio]{media/regrid_coarse.png}
+ \caption{\raggedright Mesh with \lstinline|coarsening_top_layer=True|. The top layer is remeshed with \lstinline|actual_vertical_interval_top_layer=1| vertical hex sheet.}
+ \label{fig:regrid_coarse}
+ \end{subfigure}
+ \caption{Regular mesh with one tripling layer showing the effect of \lstinline|coarsening_top_layer|.}
+ \label{fig:regrid}
\end{figure}
-\subsection{TODO: MESH for SOUTHERN CALIFORNIA}
+\subsection{TODO: Mesh for Southern California}
Mesh with high number of hexes and high resolution topography. Mesh in Parallel.
TO DO
-\subsection{TODO MESH for SANTA MONICA OVERTRUST}
+\subsection{TODO: Mesh for Santa Monica Overtrust}
Using part of the mesh for Southern California, stitching together several slices.
TO DO
-\subsection{TODO Grouping the hexes with different period resolved}
+\subsection{TODO: Grouping Hexes with Different Resolved Periods}
Checking the stability
-For execution type: \texttt{GEOCUBIT.py --meshfiles=[filename] --stability (--tomofile=[tomographic file]/--vp=[vp text value] --vs=[vs text value])}.
+For execution type: \lstinline|GEOCUBIT.py --meshfiles=[filename] --stability (--tomofile=[tomographic file]/--vp=[vp text value] --vs=[vs text value])|.
-or
+Or in the script tab of the \cubit\ GUI, type:
-in the script tab of the CUBIT GUI, type:
+\begin{lstlisting}
+>from geocubitlib.hex_metric import SEM_stability_3D # load the @\geocubit@ modules
+>mesh=SEM_stability_3D()
+>mesh.check_simulation_parameter(vp_static=vp,vs_static=vs) # test the mesh against some homogeneous velocity model
+# or
+>mesh.check_simulation_parameter(tomofile=tomofile) # test the mesh against a tomographic file
+mesh.group_timestep() # grouping hex with similar timestep
+mesh.group_period() # grouping hex with similar period
+\end{lstlisting}
-\begin{lyxcode}
->from geocubitlib.hex_metric import SEM_stability_3D $\to$ \textrm{\small{\textit{load the geocubit modules}}}\\
->mesh=SEM_stability_3D()\\
->mesh.check_simulation_parameter(vp_static=vp,vs_static=vs) $\to$ \textrm{\small{\textit{text the mesh against some homogeneous velocity model}}}\\
-or\\
->mesh.check_simulation_parameter(tomofile=tomofile) $\to$ \textrm{\small{\textit{text the mesh against a tomographic file}}}\\
-mesh.group_timestep() $\to$ \textrm{\small{\textit{grouping hex with similar timestep}}}\\
-mesh.group_period() $\to$ \textrm{\small{\textit{grouping hex with similar period}}}\\
-\end{lyxcode}
-
-
-
TO DO: picture
-
+\section{TODO: \specfem\ Mesh Format}
-
-
-\section{TODO SPECFEM3D MESH FORMAT}
-
TO DO
-\section{TODO NOTES ABOUT THE DIMENSION OF THE HEXES, THE MINIMUM PERIOD RESOLVED AND THE TIME STEP OF THE SIMULATION}
+\section{TODO: Notes about the Dimension of the Hexes, the Minimum Period Resolved and the Time Step of the Simulation}
TO DO
\bibliography{bibliography.bib}
-\end{document}
\ No newline at end of file
+\end{document}
+
Deleted: seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/example/data.zip
===================================================================
(Binary files differ)
Added: seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/example/volume_from_irregularsurf.cfg
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/example/volume_from_irregularsurf.cfg (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/example/volume_from_irregularsurf.cfg 2013-09-26 00:57:45 UTC (rev 22850)
@@ -0,0 +1,42 @@
+[cubit.options]
+cubit_info=off
+echo_info=off
+jou_info=off
+jer_info=off
+working_dir=tmp
+output_dir=output
+save_geometry_cubit = True
+save_surface_cubit = False
+
+[simulation.cpu_parameters]
+nodes = 1
+#
+[geometry.volumes]
+volume_type = layercake_volume_ascii_regulargrid_regularmap
+latitude_min = 13.90
+latitude_max = 14.20
+longitude_min = 40.65
+longitude_max = 40.95
+nx = 5
+ny = 5
+unit = geo
+# geo or utm
+
+[geometry.volumes.layercake]
+nz = 2
+#included the bottom
+bottomflat = True
+depth_bottom = -7000
+irregulargridded_surf=True
+filename = example/topoutm.dat
+geometry_format=ascii
+
+[meshing]
+map_meshing_type=regularmap
+iv_interval=5,
+size=2000
+or_mesh_scheme=map
+ntripl=1
+smoothing=False
+coarsening_top_layer=False
+tripl=2,
\ No newline at end of file
Modified: seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/boundary_definition.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/boundary_definition.py 2013-09-26 00:49:25 UTC (rev 22849)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/boundary_definition.py 2013-09-26 00:57:45 UTC (rev 22850)
@@ -1,3 +1,26 @@
+#############################################################################
+# boundary_definition.py
+# this file is part of GEOCUBIT #
+# #
+# Created by Emanuele Casarotti #
+# Copyright (c) 2008 Istituto Nazionale di Geofisica e Vulcanologia #
+# #
+#############################################################################
+# #
+# GEOCUBIT is free software: you can redistribute it and/or modify #
+# it under the terms of the GNU General Public License as published by #
+# the Free Software Foundation, either version 3 of the License, or #
+# (at your option) any later version. #
+# #
+# GEOCUBIT is distributed in the hope that it will be useful, #
+# but WITHOUT ANY WARRANTY; without even the implied warranty of #
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
+# GNU General Public License for more details. #
+# #
+# You should have received a copy of the GNU General Public License #
+# along with GEOCUBIT. If not, see <http://www.gnu.org/licenses/>. #
+# #
+#############################################################################
try:
import start as start
cubit = start.start_cubit()
@@ -8,9 +31,7 @@
print 'error importing cubit, check if cubit is installed'
pass
-def list2str(l):
- if not isinstance(l,list): l=list(l)
- return ' '.join(str(x) for x in l)
+from utilities import list2str
def map_boundary(cpuxmin=0,cpuxmax=1,cpuymin=0,cpuymax=1,cpux=1,cpuy=1):
@@ -95,13 +116,21 @@
def define_surf(ip=0,cpuxmin=0,cpuxmax=1,cpuymin=0,cpuymax=1,cpux=1,cpuy=1):
"""
- define the surfaces defining the boundaries of the volume
+ define the absorbing surfaces for a layered topological box where boundary are surfaces parallel to the axis.
+ it returns absorbing_surf,absorbing_surf_xmin,absorbing_surf_xmax,absorbing_surf_ymin,absorbing_surf_ymax,absorbing_surf_bottom,topo_surf
+ where
+ absorbing_surf is the list of all the absorbing boundary surf
+ absorbing_surf_xmin is the list of the absorbing boundary surfaces that correnspond to x=xmin
+ ...
+ absorbing_surf_bottom is the list of the absorbing boundary surfaces that correspond to z=zmin
"""
+ from utilities import get_v_h_list
#
from sets import Set
def product(*args, **kwds):
# product('ABCD', 'xy') --> Ax Ay Bx By Cx Cy Dx Dy
# product(range(2), repeat=3) --> 000 001 010 011 100 101 110 111
+ # for compatibility with python2.5
pools = map(tuple, args) * kwds.get('repeat', 1)
result = [[]]
for pool in pools:
@@ -123,6 +152,11 @@
ymin_box=cubit.get_total_bounding_box("volume",list_vol)[3]
ymax_box=cubit.get_total_bounding_box("volume",list_vol)[4]
list_surf=cubit.parse_cubit_list("surface","all")
+
+ absorbing_surface_distance_tolerance=0.001
+ topographic_surface_distance_tolerance=0.1
+ topographic_surface_normal_tolerance=0.4
+
lv=[]
for k in list_surf:
sbox=cubit.get_bounding_box('surface',k)
@@ -140,16 +174,19 @@
dzmin=abs(sbox[6] - zmin_box)/max(abs(sbox[6]),abs(zmin_box))
normal=cubit.get_surface_normal(k)
zn=normal[2]
- if dzmax <= 0.1 and zn > 0.4:
+ if dzmax <= topographic_surface_distance_tolerance and zn > topographic_surface_normal_tolerance:
top_surf.append(k)
list_vertex=cubit.get_relatives('surface',k,'vertex')
for v in list_vertex:
valence=cubit.get_valence(v)
if valence <= 4: #valence 3 is a corner, 4 is a vertex between 2 volumes, > 4 is a vertex not in the boundaries
lv.append(v)
- elif dzmin <= 0.001 and zn < -0.7:
+ elif dzmin <= 0.001 and zn < -1+topographic_surface_normal_tolerance:
bottom_surf.append(k)
+ if len(top_surf) ==0: #assuming that one topo surface need to be selected
+ _,_,_,_,_,top_surf=get_v_h_list(list_vol,chktop=False)
lp=[]
+ labelp=[]
combs=product(lv,lv)
for comb in combs:
v1=comb[0]
@@ -157,19 +194,36 @@
c=Set(cubit.get_relatives("vertex",v1,"curve")) & Set(cubit.get_relatives("vertex",v2,"curve"))
if len(c) == 1:
p=cubit.get_center_point("curve",list(c)[0])
- lp.append(p)
+ labelp.append(list(c)[0])
+ labelps=Set(labelp)
+ for c in labelps:
+ p=cubit.get_center_point("curve",c)
+ lp.append(p)
+
for k in list_surf:
center_point = cubit.get_center_point("surface", k)
for p in lp:
- if abs((center_point[0] - p[0])/p[0]) <= 0.001 and abs((center_point[1] - p[1])/p[1]) <= 0.001:
- absorbing_surf.append(k)
- break
+ try:
+ if abs((center_point[0] - p[0])/p[0]) <= absorbing_surface_distance_tolerance and abs((center_point[1] - p[1])/p[1]) <= absorbing_surface_distance_tolerance:
+ absorbing_surf.append(k)
+ break
+ except:
+ if -1 <= center_point[0] <= 1 and -1 <= center_point[1] <= 1:
+ absorbing_surf.append(k)
+ break
#
four_side=True
if four_side:
- xmin,ymin,xmax,ymax=define_4side_lateral_surfaces()
- abs_xmin,abs_xmax,abs_ymin,abs_ymax=lateral_boundary_are_absorbing(ip,cpuxmin,cpuxmax,cpuymin,cpuymax,cpux,cpuy)
-
+ xmintmp,ymintmp,xmaxtmp,ymaxtmp=define_4side_lateral_surfaces()
+ xmin=list(Set(xmintmp)-Set(xmaxtmp))
+ xmax=list(Set(xmaxtmp)-Set(xmintmp))
+ ymin=list(Set(ymintmp)-Set(ymaxtmp))
+ ymax=list(Set(ymaxtmp)-Set(ymintmp))
+ abs_xmintmp,abs_xmaxtmp,abs_ymintmp,abs_ymaxtmp=lateral_boundary_are_absorbing(ip,cpuxmin,cpuxmax,cpuymin,cpuymax,cpux,cpuy)
+ abs_xmin=list(Set(abs_xmintmp)-Set(abs_xmaxtmp))
+ abs_xmax=list(Set(abs_xmaxtmp)-Set(abs_xmintmp))
+ abs_ymin=list(Set(abs_ymintmp)-Set(abs_ymaxtmp))
+ abs_ymax=list(Set(abs_ymaxtmp)-Set(abs_ymintmp))
return absorbing_surf,abs_xmin,abs_xmax,abs_ymin,abs_ymax,top_surf,bottom_surf,xmin,ymin,xmax,ymax
@@ -180,9 +234,20 @@
list_name=map(lambda x: 'vol'+x,map(str,list_vol))
return list_vol,list_name
-def build_block(vol_list,name,id_0=1):
+def build_block(vol_list,name,id_0=1,top_surf=None,optionsea=False):
#
from sets import Set
+ if optionsea:
+ sea=optionsea['sea']
+ seaup=optionsea['seaup']
+ sealevel=optionsea['sealevel']
+ seathres=optionsea['seathres']
+ else:
+ sea=False
+ seaup=False
+ sealevel=False
+ seathres=False
+
#
block_list=cubit.get_block_id_list()
if len(block_list) > 0:
@@ -193,13 +258,35 @@
id_block+=1
v_other=Set(vol_list)-Set([v])
#command= 'block '+str(id_block)+' hex in node in vol '+str(v)+' except hex in vol '+str(list(v_other))
- command= 'block '+str(id_block)+' hex in vol '+str(v)+' except hex in vol '+str(list(v_other))
- print command
- command = command.replace("["," ").replace("]"," ")
- cubit.cmd(command)
- command = "block "+str(id_block)+" name '"+n+"'"
- cubit.cmd(command)
+ if sea and v == vol_list[-1]:
+ cubit.cmd('set duplicate block elements off')
+ tsurf_string=" ".join(str(x) for x in top_surf)
+ #sea
+ command= 'block '+str(id_block)+' hex in node in surf '+tsurf_string+' with Z_coord < '+str(seathres)
+ cubit.cmd(command)
+ command = "block "+str(id_block)+" name 'sea"+n+"'"
+ cubit.cmd(command)
+ if not seaup:
+ id_block+=1
+ command= 'block '+str(id_block)+' hex in node in surf '+tsurf_string+' with (Z_coord > '+str(seathres)+' and Z_coord < '+str(sealevel)
+ cubit.cmd(command)
+ command = "block "+str(id_block)+" name 'shwater"+n+"'"
+ cubit.cmd(command)
+ id_block+=1
+ command= 'block '+str(id_block)+' hex in node in surf '+tsurf_string+' with Z_coord >= '+str(sealevel)
+ cubit.cmd(command)
+ command = "block "+str(id_block)+" name 'continent"+n+"'"
+ cubit.cmd(command)
+ else:
+ command= 'block '+str(id_block)+' hex in vol '+str(v)+' except hex in vol '+str(list(v_other))
+ print command
+ command = command.replace("["," ").replace("]"," ")
+ cubit.cmd(command)
+ command = "block "+str(id_block)+" name '"+n+"'"
+ cubit.cmd(command)
+
+
def build_block_side(surf_list,name,obj='surface',id_0=1):
id_nodeset=cubit.get_next_nodeset_id()
id_block=id_0
@@ -246,26 +333,13 @@
cpuy=keys.get("cpuy",1)
cpuxmax=keys.get("cpuxmax",cpux)
cpuymax=keys.get("cpuymax",cpuy)
+ optionsea=keys.get("optionsea",False)
#
if parallel:
absorbing_surf,abs_xmin,abs_xmax,abs_ymin,abs_ymax,top_surf,bottom_surf,xmin,ymin,xmax,ymax=define_surf(ip=ip,cpuxmin=cpuxmin,cpuxmax=cpuxmax,cpuymin=cpuymin,cpuymax=cpuymax,cpux=cpux,cpuy=cpuy)
- #
- #id_side=cubit.get_next_block_id()
- #try:
- # entities=args[0]
- #except:
- # entities=['face']
- #for entity in entities:
- # build_block_side(top_surf,entity+'_topo',obj=entity,id_0=1) #topo is block 1 so id_0=1
- # build_block_side(bottom_surf,entity+'_abs_bottom',obj=entity,id_0=2)
- # build_block_side(xmin,entity+'_xmin',obj=entity,id_0=3)
- # build_block_side(ymin,entity+'_ymin',obj=entity,id_0=4)
- # build_block_side(xmax,entity+'_xmax',obj=entity,id_0=5)
- # build_block_side(ymax,entity+'_ymax',obj=entity,id_0=6)
- #
id_0=cubit.get_next_block_id()
v_list,name_list=define_block()
- build_block(v_list,name_list,id_0)
+ build_block(v_list,name_list,id_0,top_surf,optionsea=optionsea)
#
elif closed:
surf=define_absorbing_surf_sphere()
@@ -331,14 +405,24 @@
icurvestr=str(icurve)
orient_nodes_surf=[]
#
- cubit.cmd('del group sl')
+ #get the nodes on a surface, I don't use the method get_surface_nodes since it has different behavior in cubit12.2 and cubit13.2+
+ k=cubit.get_id_from_name('sl')
+ if k!=0:
+ cubit.cmd('del group sl')
+ else:
+ print 'initializing group sl'
cubit.cmd("group 'sl' add node in surf "+lsurf)
group1 = cubit.get_id_from_name("sl")
nodes_ls =list(cubit.get_group_nodes(group1))
nnode=len(nodes_ls)
#
+ #get the nodes on curves
orient=[]
- cubit.cmd('del group n1')
+ k=cubit.get_id_from_name('n1')
+ if k!=0:
+ cubit.cmd('del group n1')
+ else:
+ print 'initializing group n1'
cubit.cmd("group 'n1' add node in curve "+icurvestr)
x=cubit.get_bounding_box('curve', icurve)
if x[2]>x[5]:
@@ -392,9 +476,13 @@
if length[-1][1] == 3:
curve_vertical.append(l)
#
- icurve=list2str(curve_vertical)
- cubit.cmd('del group curve_vertical')
- cubit.cmd("group 'curve_vertical' add node in curve "+icurve)
+ kcurve=list2str(curve_vertical)
+ k=cubit.get_id_from_name('curve_vertical')
+ if k!=0:
+ cubit.cmd('del group curve_vertical')
+ else:
+ print 'initializing group curve_vertical'
+ cubit.cmd("group 'curve_vertical' add node in curve "+kcurve)
group1 = cubit.get_id_from_name('curve_vertical')
nodes_curve = list(cubit.get_group_nodes(group1))
for n in nodes_curve:
@@ -491,9 +579,11 @@
boundary['node_curve_xminymax']=c_xminymax
boundary['node_curve_xmaxymin']=c_xmaxymin
boundary['node_curve_xmaxymax']=c_xmaxymax
- #
- #
- #
+
+ #print boundary['node_curve_xminymin']
+ #print boundary['node_curve_xminymax']
+ #print boundary['node_curve_xmaxymin']
+ #print boundary['node_curve_xmaxymax']
entities=['face']
#
for entity in entities:
@@ -517,14 +607,24 @@
block=3 #change here..... must be 1 @@@@@@@@@
ty=None
ty=cubit.get_block_element_type(block)
- if ty != 'HEX8': cubit.cmd('del block '+str(block))
+ if ty=='':
+ pass
+ elif ty == 'HEX8':
+ pass
+ else:
+ cubit.cmd('del block '+str(block))
build_block_side(top_surf,refname,obj=entity,id_0=1001)
#
refname=entity+'_bottom'
block=4 #change here..... must be 2 @@@@@@@@
ty=None
ty=cubit.get_block_element_type(block)
- if ty != 'HEX8': cubit.cmd('del block '+str(block))
+ if ty=='':
+ pass
+ elif ty == 'HEX8':
+ pass
+ else:
+ cubit.cmd('del block '+str(block))
build_block_side(bottom_surf,refname,obj=entity,id_0=1002)
#
#
Modified: seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/cubit2specfem3d.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/cubit2specfem3d.py 2013-09-26 00:49:25 UTC (rev 22849)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/cubit2specfem3d.py 2013-09-26 00:57:45 UTC (rev 22850)
@@ -126,6 +126,8 @@
print 'error importing cubit, check if cubit is installed'
pass
+from utilities import get_cubit_version
+
class mtools(object):
def __init__(self,frequency,list_surf,list_vp):
super(mtools, self).__init__()
@@ -340,28 +342,47 @@
return self.ddt[0],self.dr[0]
class mesh(object,mesh_tools):
- def __init__(self):
+ def __init__(self,hex27=False,cpml=False,cpml_size=False,top_absorbing=False):
super(mesh, self).__init__()
self.mesh_name='mesh_file'
self.nodecoord_name='nodes_coords_file'
self.material_name='materials_file'
self.nummaterial_name='nummaterial_velocity_file'
self.absname='absorbing_surface_file'
- self.freename='free_surface_file'
+ self.cpmlname='absorbing_cpml_file'
+ self.freename='free_or_absorbing_surface_file_zmax'
self.recname='STATIONS'
- version_cubit=float(cubit.get_version())
+ version_cubit=get_cubit_version()
if version_cubit >= 12:
self.face='SHELL4'
else:
self.face='QUAD4'
self.hex='HEX'
+ self.hex27=hex27
self.edge='BAR2'
self.topo='face_topo'
+ self.topography=None
+ self.free=None
+ self.freetxt='free'
self.rec='receivers'
+ self.cpml=cpml
+ if cpml:
+ if cpml_size:
+ self.size=cpml_size
+ else:
+ print 'please specify cmpl size if you want to use cpml'
+ self.top_absorbing=top_absorbing
+ if hex27: cubit.cmd('block all except block 1001 1002 1003 1004 1005 1006 element type hex27')
self.block_definition()
self.ngll=5
self.percent_gll=0.172
self.point_wavelength=5
+ self.xmin=False
+ self.ymin=False
+ self.zmin=False
+ self.xmax=False
+ self.ymax=False
+ self.zmax=False
cubit.cmd('compress all')
def __repr__(self):
pass
@@ -394,16 +415,13 @@
imaterial = 3
else :
imaterial = 0
-
-## DK DK added this default value for all parameters to avoid undefined parameters for older meshes in the "examples" directory that do not contain Qkappa
- vel,vs,rho,qk,qmu,ani=(0,0,0,9999.,9999.,0)
-
+ #
if nattrib > 1:
# material flag:
# positive => material properties,
# negative => interface/tomography domain
flag=int(cubit.get_block_attribute_value(block,0))
- if flag > 0 and nattrib >= 2:
+ if 0< flag and nattrib >= 2:
vel=cubit.get_block_attribute_value(block,1)
if nattrib >= 3:
vs=cubit.get_block_attribute_value(block,2)
@@ -430,7 +448,7 @@
kind='tomography'
elif nattrib == 1:
flag=cubit.get_block_attribute_value(block,0)
- print 'only 1 attribute ', name,block,flag
+ #print 'only 1 attribute ', name,block,flag
vel,vs,rho,qk,qmu,ani=(0,0,0,9999.,9999.,0)
else:
flag=block
@@ -451,7 +469,10 @@
block_bc_flag.append(4)
block_bc.append(block)
bc[block]=4 #face has connectivity = 4
- if name == self.topo or block == 1001: topography_face=block
+ if name == self.topo or block == 1001:
+ self.topography=block
+ if self.freetxt in name:
+ self.free=block
elif ty == 'SPHERE':
pass
else:
@@ -479,13 +500,7 @@
else:
print 'nodeset '+name+' not defined'
self.receivers=None
- print block_mat
- print block_flag
- print block_bc
- print block_bc_flag
- print material
- print bc
- print topography_face
+
try:
self.block_mat=block_mat
self.block_flag=block_flag
@@ -493,7 +508,16 @@
self.block_bc_flag=block_bc_flag
self.material=material
self.bc=bc
- self.topography=topography_face
+ print 'HEX Blocks:'
+ for m,f in zip(self.block_mat,self.block_flag):
+ print 'block ',m,'material flag ',f
+ print 'Absorbing Boundary Conditions:'
+ for m,f in zip(self.block_bc,self.block_bc_flag):
+ print 'bc ',m,'bc flag ',f
+ print 'Topography (free surface)'
+ print self.topography
+ print 'Free surface'
+ print self.free
except:
print '****************************************'
print 'sorry, no blocks or blocks not properly defined'
@@ -503,14 +527,40 @@
print block_bc_flag
print material
print bc
- print topography
print '****************************************'
+
+
+
+
+
+
+ def get_hex_connectivity(self,ind):
+ if self.hex27:
+ cubit.silent_cmd('group "nh" add Node in hex '+str(ind))
+ group1 = cubit.get_id_from_name("nh")
+ result=cubit.get_group_nodes(group1)
+ cubit.cmd('del group '+str(group1))
+ else:
+ result=cubit.get_connectivity('hex',ind)
+ return result
+ #
+ def get_face_connectivity(self,ind):
+ if self.hex27:
+ cubit.silent_cmd('group "nf" add Node in face '+str(ind))
+ group1 = cubit.get_id_from_name("nf")
+ result=cubit.get_group_nodes(group1)
+ cubit.cmd('del group '+str(group1))
+ else:
+ result=cubit.get_connectivity('face',ind)
+ return result
+
+
def mat_parameter(self,properties):
- print properties
+ #print properties
#format nummaterials file: #material_domain_id #material_id #rho #vp #vs #Q_kappa #Q_mu #anisotropy_flag
imaterial=properties[0]
flag=properties[1]
- print 'prop',flag
+ #print 'prop',flag
if flag > 0:
vel=properties[2]
if properties[2] is None and type(vel) != str:
@@ -540,7 +590,7 @@
helpstring="#material_domain_id #material_id #rho #vp #vs #Q_kappa #Q_mu #anisotropy"
txt='%1i %3i %s \n' % (properties[0],properties[1],helpstring)
else:
- helpstring=" --> syntax: #material_domain_id #material_id #rho #vp #vs #Q_kappa #Q_mu #anisotropy"
+ helpstring=" --> sintax: #material_domain_id #material_id #rho #vp #vs #Q_kappa #Q_mu #anisotropy"
txt='%1i %3i %s %s\n' % (properties[0],properties[1],properties[2],helpstring)
elif flag < 0:
if properties[2] == 'tomography':
@@ -548,10 +598,10 @@
elif properties[2] == 'interface':
txt='%1i %3i %s %s %1i %1i\n' % (properties[0],properties[1],properties[2],properties[3],properties[4],properties[5])
else:
- helpstring=" --> syntax: #material_domain_id 'tomography' #file_name "
+ helpstring=" --> sintax: #material_domain_id 'tomography' #file_name "
txt='%1i %3i %s %s \n' % (properties[0],properties[1],properties[2],helpstring)
#
- print txt
+ #print txt
return txt
def nummaterial_write(self,nummaterial_name):
print 'Writing '+nummaterial_name+'.....'
@@ -560,32 +610,78 @@
#name=cubit.get_exodus_entity_name('block',block)
nummaterial.write(str(self.mat_parameter(self.material[block])))
nummaterial.close()
+
+ def create_hexnode_string(self,hexa):
+ nodes=self.get_hex_connectivity(hexa)
+ #nodes=self.jac_check(nodes) #is it valid for 3D? TODO
+ if self.hex27:
+ ordered_nodes=[hexa]+list(nodes[:20])+[nodes[21]]+[nodes[25]]+[nodes[24]]+[nodes[26]]+[nodes[23]]+[nodes[22]]+[nodes[20]]
+ txt=' '.join(str(x) for x in ordered_nodes)
+ txt=txt+'\n'
+ #txt=('%10i %10i %10i %10i %10i %10i %10i %10i ')% nodes[:8] #first 8 nodes following specfem3d numbering convenction..
+ #txt=txt+('%10i %10i %10i %10i %10i %10i %10i %10i ')% nodes[8:16] #middle 12 nodes following specfem3d numbering convenction..
+ #txt=txt+('%10i %10i %10i %10i ')% nodes[16:20]
+ #txt=txt+('%10i %10i %10i %10i %10i %10i ')% (nodes[21], nodes[25], nodes[24], nodes[26], nodes[23], nodes[22])
+ #txt=txt+('%10i\n ')% nodes[20] #center volume
+ else:
+ txt=str(hexa)+' '+' '.join(str(x) for x in nodes)
+ txt=txt+'\n'
+ #txt=('%10i %10i %10i %10i %10i %10i %10i %10i\n')% nodes[:]
+ return txt
+
+ def create_facenode_string(self,hexa,face,normal=None,cknormal=True):
+ nodes=self.get_face_connectivity(face)
+ if cknormal:
+ nodes_ok=self.normal_check(nodes[0:4],normal)
+ if self.hex27: nodes_ok2=self.normal_check(nodes[4:8],normal)
+ else:
+ nodes_ok=nodes[0:4]
+ if self.hex27: nodes_ok2=nodes[4:8]
+ #
+ if self.hex27:
+ ordered_nodes=[hexa]+list(nodes_ok)+list(nodes_ok2)+[nodes[8]]
+ txt=' '.join(str(x) for x in ordered_nodes)
+ txt=txt+'\n'
+ #txt=('%10i %10i %10i %10i %10i ') % (hexa,nodes_ok[0],nodes_ok[1],nodes_ok[2],nodes_ok[3]) #first 4 nodes following specfem3d numbering convenction..
+ #txt=txt+('%10i %10i %10i %10i ')% (nodes_ok2[0],nodes_ok2[1],nodes_ok2[2],nodes_ok2[3]) #middle 4 nodes following specfem3d numbering convenction..
+ #txt=txt+('%10i\n')% nodes[8]
+ else:
+ txt=str(hexa)+' '+' '.join(str(x) for x in nodes_ok)
+ txt=txt+'\n'
+ #txt=('%10i %10i %10i %10i %10i\n') % (hexa,nodes_ok[0],nodes_ok[1],nodes_ok[2],nodes_ok[3])
+ return txt
+
+
def mesh_write(self,mesh_name):
meshfile=open(mesh_name,'w')
print 'Writing '+mesh_name+'.....'
num_elems=cubit.get_hex_count()
- print ' number of elements:',str(num_elems)
+ print ' total number of elements:',str(num_elems)
meshfile.write(str(num_elems)+'\n')
for block,flag in zip(self.block_mat,self.block_flag):
- #print block,flag
hexes=cubit.get_block_hexes(block)
- #print len(hexes)
+ print 'block ',block,' hexes ',len(hexes)
for hexa in hexes:
- #print hexa
- nodes=cubit.get_connectivity('Hex',hexa)
- #nodes=self.jac_check(nodes) #is it valid for 3D? TODO
- txt=('%10i ')% hexa
- txt=txt+('%10i %10i %10i %10i %10i %10i %10i %10i\n')% nodes[:]
+ txt=self.create_hexnode_string(hexa)
meshfile.write(txt)
meshfile.close()
def material_write(self,mat_name):
mat=open(mat_name,'w')
print 'Writing '+mat_name+'.....'
for block,flag in zip(self.block_mat,self.block_flag):
+ print 'block ',block,'flag ',flag
hexes=cubit.get_block_hexes(block)
for hexa in hexes:
mat.write(('%10i %10i\n') % (hexa,flag))
mat.close()
+ def get_extreme(self,c,cmin,cmax):
+ if not cmin and not cmax:
+ cmin=c
+ cmax=c
+ else:
+ if c<cmin: cmin=c
+ if c>cmax: cmax=c
+ return cmin,cmax
def nodescoord_write(self,nodecoord_name):
nodecoord=open(nodecoord_name,'w')
print 'Writing '+nodecoord_name+'.....'
@@ -597,6 +693,9 @@
for node in node_list:
x,y,z=cubit.get_nodal_coordinates(node)
+ self.xmin,self.xmax=self.get_extreme(x,self.xmin,self.xmax)
+ self.ymin,self.ymax=self.get_extreme(y,self.ymin,self.ymax)
+ self.zmin,self.zmax=self.get_extreme(z,self.zmin,self.zmax)
txt=('%10i %20f %20f %20f\n') % (node,x,y,z)
nodecoord.write(txt)
nodecoord.close()
@@ -614,7 +713,7 @@
for block,flag in zip(self.block_bc,self.block_bc_flag):
if block == self.topography:
name=cubit.get_exodus_entity_name('block',block)
- print ' block name:',name,'id:',block
+ print 'free surface (topography) block name:',name,'id:',block
quads_all=cubit.get_block_faces(block)
print ' number of faces = ',len(quads_all)
dic_quads_all=dict(zip(quads_all,quads_all))
@@ -625,13 +724,114 @@
for f in faces:
if dic_quads_all.has_key(f):
#print f
- nodes=cubit.get_connectivity('face',f)
- nodes_ok=self.normal_check(nodes,normal)
- txt='%10i %10i %10i %10i %10i\n' % (h,nodes_ok[0],nodes_ok[1],nodes_ok[2],nodes_ok[3])
+ txt=self.create_facenode_string(h,f,normal,cknormal=True)
freehex.write(txt)
- freehex.close()
+ freehex.close()
+ elif block == self.free:
+ name=cubit.get_exodus_entity_name('block',block)
+ print 'free surface block name:',name,'id:',block
+ quads_all=cubit.get_block_faces(block)
+ print ' number of faces = ',len(quads_all)
+ dic_quads_all=dict(zip(quads_all,quads_all))
+ freehex.write('%10i\n' % len(quads_all))
+ list_hex=cubit.parse_cubit_list('hex','all')
+ for h in list_hex:
+ faces=cubit.get_sub_elements('hex',h,2)
+ for f in faces:
+ if dic_quads_all.has_key(f):
+ txt=self.create_facenode_string(h,f,normal,cknormal=False)
+ freehex.write(txt)
+ freehex.close()
cubit.cmd('set info on')
cubit.cmd('set echo on')
+ def check_cmpl_size(self,case='x'):
+ if case=='x':
+ vmaxtmp=self.xmax
+ vmintmp=self.xmin
+ elif case=='y':
+ vmaxtmp=self.ymax
+ vmintmp=self.ymin
+ elif case=='z':
+ vmaxtmp=self.zmax
+ vmintmp=self.zmin
+
+ if self.size > .3*(vmaxtmp-vmintmp):
+ print 'please select the size of cpml less than 30% of the '+case+' size of the volume'
+ print vmaxtmp-vmintmp,.3*(vmaxtmp-vmintmp)
+ print 'cmpl set to false, no '+self.cpmlname+' file will be created'
+ return False,False
+ else:
+ vmin=vmintmp+self.size
+ vmax=vmaxtmp-self.size
+ return vmin,vmax
+ def select_cpml(self):
+ xmin,xmax=self.check_cmpl_size(case='x')
+ ymin,ymax=self.check_cmpl_size(case='y')
+ zmin,zmax=self.check_cmpl_size(case='z')
+ #
+ if xmin is False or xmax is False or ymin is False or ymax is False or zmin is False or zmax is False:
+ return False
+ else:
+ txt="group 'hxmin' add hex with X_coord < "+str(xmin)
+ cubit.cmd(txt)
+ txt="group 'hxmax' add hex with X_coord > "+str(xmax)
+ cubit.cmd(txt)
+ txt="group 'hymin' add hex with Y_coord < "+str(ymin)
+ cubit.cmd(txt)
+ txt="group 'hymax' add hex with Y_coord > "+str(ymax)
+ cubit.cmd(txt)
+ txt="group 'hzmin' add hex with Z_coord < "+str(zmin)
+ cubit.cmd(txt)
+ txt="group 'hzmax' add hex with Z_coord > "+str(zmax)
+ cubit.cmd(txt)
+ from sets import Set
+ group1 = cubit.get_id_from_name("hxmin")
+ cpml_xmin =Set(list(cubit.get_group_hexes(group1)))
+ group1 = cubit.get_id_from_name("hymin")
+ cpml_ymin =Set(list(cubit.get_group_hexes(group1)))
+ group1 = cubit.get_id_from_name("hxmax")
+ cpml_xmax =Set(list(cubit.get_group_hexes(group1)))
+ group1 = cubit.get_id_from_name("hymax")
+ cpml_ymax =Set(list(cubit.get_group_hexes(group1)))
+ group1 = cubit.get_id_from_name("hzmin")
+ cpml_zmin =Set(list(cubit.get_group_hexes(group1)))
+ if self.top_absorbing:
+ group1 = cubit.get_id_from_name("hzmax")
+ cpml_zmax =Set(list(cubit.get_group_hexes(group1)))
+ else:
+ cpml_zmax =Set([])
+ cpml_all=cpml_ymin | cpml_ymax | cpml_xmin | cpml_xmax | cpml_zmin | cpml_zmax
+ cpml_x=cpml_all-cpml_zmin-cpml_ymin-cpml_ymax-cpml_zmax
+ cpml_y=cpml_all-cpml_zmin-cpml_xmin-cpml_xmax-cpml_zmax
+ cpml_xy=cpml_all-cpml_zmin-cpml_y-cpml_x-cpml_zmax
+ cpml_z=cpml_all-cpml_xmin-cpml_ymin-cpml_ymax-cpml_xmax
+ cpml_xz=cpml_zmin-cpml_ymin-cpml_ymax-cpml_z
+ cpml_yz=cpml_zmin-cpml_xmin-cpml_xmax-cpml_z
+ cpml_xyz=cpml_zmin-cpml_xz-cpml_yz-cpml_z
+ txt=' '.join(str(h) for h in cpml_x)
+ cubit.cmd("group 'x_cpml' add hex "+txt)
+ txt=' '.join(str(h) for h in cpml_y)
+ cubit.cmd("group 'y_cpml' add hex "+txt)
+ txt=' '.join(str(h) for h in cpml_z)
+ cubit.cmd("group 'z_cpml' add hex "+txt)
+ txt=' '.join(str(h) for h in cpml_xy)
+ cubit.cmd("group 'xy_cpml' add hex "+txt)
+ txt=' '.join(str(h) for h in cpml_xz)
+ cubit.cmd("group 'xz_cpml' add hex "+txt)
+ txt=' '.join(str(h) for h in cpml_yz)
+ cubit.cmd("group 'yz_cpml' add hex "+txt)
+ txt=' '.join(str(h) for h in cpml_xyz)
+ cubit.cmd("group 'xyz_cpml' add hex "+txt)
+ return cpml_x,cpml_y,cpml_z,cpml_xy,cpml_xz,cpml_yz,cpml_xyz
+
+
+
+
+
+
+
+
+
def abs_write(self,absname=None):
import re
cubit.cmd('set info off')
@@ -639,86 +839,108 @@
cubit.cmd('set journal off')
from sets import Set
if not absname: absname=self.absname
- #
- #
- list_hex=cubit.parse_cubit_list('hex','all')
- for block,flag in zip(self.block_bc,self.block_bc_flag):
- if block != self.topography:
- name=cubit.get_exodus_entity_name('block',block)
- print ' block name:',name,'id:',block
- cknormal=True
- if re.search('xmin',name):
- print 'xmin'
- abshex_local=open(absname+'_xmin','w')
- normal=(-1,0,0)
- elif re.search('xmax',name):
- print "xmax"
- abshex_local=open(absname+'_xmax','w')
- normal=(1,0,0)
- elif re.search('ymin',name):
- print "ymin"
- abshex_local=open(absname+'_ymin','w')
- normal=(0,-1,0)
- elif re.search('ymax',name):
- print "ymax"
- abshex_local=open(absname+'_ymax','w')
- normal=(0,1,0)
- elif re.search('bottom',name):
- print "bottom"
- abshex_local=open(absname+'_bottom','w')
- normal=(0,0,-1)
- elif re.search('abs',name):
- print "abs all - no implemented yet"
- cknormal=False
- abshex_local=open(absname,'w')
- else:
- if block == 1003:
+
+ if self.cpml:
+ if not absname: absname=self.cpmlname
+ print 'Writing cpml'+absname+'.....'
+ list_cpml=self.select_cpml()
+ if list_cpml is False:
+ print 'error writing cpml files'
+ return
+ else:
+ abshex_cpml=open(absname,'w')
+ hexcount=sum(map(len,list_cpml))
+ abshex_cpml.write(('%10i\n') % (hexcount))
+ for icpml,lcpml in enumerate(list_cpml):
+ for hexa in lcpml:
+ abshex_cpml.write(('%10i %10i\n') % (hexa,icpml))
+
+
+ stacey_absorb=True
+ if stacey_absorb:
+ #
+ #
+ if not absname: absname=self.absname
+ list_hex=cubit.parse_cubit_list('hex','all')
+ for block,flag in zip(self.block_bc,self.block_bc_flag):
+ if block != self.topography:
+ name=cubit.get_exodus_entity_name('block',block)
+ print ' block name:',name,'id:',block
+ cknormal=True
+ abshex_local=False
+ if re.search('xmin',name):
print 'xmin'
abshex_local=open(absname+'_xmin','w')
normal=(-1,0,0)
- elif block == 1004:
+ elif re.search('xmax',name):
+ print "xmax"
+ abshex_local=open(absname+'_xmax','w')
+ normal=(1,0,0)
+ elif re.search('ymin',name):
print "ymin"
abshex_local=open(absname+'_ymin','w')
normal=(0,-1,0)
- elif block == 1005:
- print "xmax"
- abshex_local=open(absname+'_xmax','w')
- normal=(1,0,0)
- elif block == 1006:
+ elif re.search('ymax',name):
print "ymax"
abshex_local=open(absname+'_ymax','w')
normal=(0,1,0)
- elif block == 1002:
+ elif re.search('bottom',name):
print "bottom"
abshex_local=open(absname+'_bottom','w')
normal=(0,0,-1)
- #
- #
- quads_all=cubit.get_block_faces(block)
- dic_quads_all=dict(zip(quads_all,quads_all))
- print ' number of faces = ',len(quads_all)
- abshex_local.write('%10i\n' % len(quads_all))
- #command = "group 'list_hex' add hex in face "+str(quads_all)
- #command = command.replace("["," ").replace("]"," ").replace("("," ").replace(")"," ")
- #cubit.cmd(command)
- #group=cubit.get_id_from_name("list_hex")
- #list_hex=cubit.get_group_hexes(group)
- #command = "delete group "+ str(group)
- #cubit.cmd(command)
- for h in list_hex:
- faces=cubit.get_sub_elements('hex',h,2)
- for f in faces:
- if dic_quads_all.has_key(f):
- nodes=cubit.get_connectivity('face',f)
- if cknormal:
- nodes_ok=self.normal_check(nodes,normal)
- else:
- nodes_ok=nodes
- txt='%10i %10i %10i %10i %10i\n' % (h,nodes_ok[0],nodes_ok[1],nodes_ok[2],nodes_ok[3])
- abshex_local.write(txt)
- abshex_local.close()
- cubit.cmd('set info on')
- cubit.cmd('set echo on')
+ elif re.search('abs',name):
+ print "abs all - experimental, check the output"
+ cknormal=False
+ abshex_local=open(absname,'w')
+ else:
+ if block == 1003:
+ print 'xmin'
+ abshex_local=open(absname+'_xmin','w')
+ normal=(-1,0,0)
+ elif block == 1004:
+ print "ymin"
+ abshex_local=open(absname+'_ymin','w')
+ normal=(0,-1,0)
+ elif block == 1005:
+ print "xmax"
+ abshex_local=open(absname+'_xmax','w')
+ normal=(1,0,0)
+ elif block == 1006:
+ print "ymax"
+ abshex_local=open(absname+'_ymax','w')
+ normal=(0,1,0)
+ elif block == 1002:
+ print "bottom"
+ abshex_local=open(absname+'_bottom','w')
+ normal=(0,0,-1)
+ elif block == 1000:
+ print "custumized"
+ abshex_local=open(absname,'w')
+ cknormal=False
+ normal=None
+ #
+ #
+ if abshex_local:
+ quads_all=cubit.get_block_faces(block)
+ dic_quads_all=dict(zip(quads_all,quads_all))
+ print ' number of faces = ',len(quads_all)
+ abshex_local.write('%10i\n' % len(quads_all))
+ #command = "group 'list_hex' add hex in face "+str(quads_all)
+ #command = command.replace("["," ").replace("]"," ").replace("("," ").replace(")"," ")
+ #cubit.cmd(command)
+ #group=cubit.get_id_from_name("list_hex")
+ #list_hex=cubit.get_group_hexes(group)
+ #command = "delete group "+ str(group)
+ #cubit.cmd(command)
+ for h in list_hex:
+ faces=cubit.get_sub_elements('hex',h,2)
+ for f in faces:
+ if dic_quads_all.has_key(f):
+ txt=self.create_facenode_string(h,f,normal=normal,cknormal=cknormal)
+ abshex_local.write(txt)
+ abshex_local.close()
+ cubit.cmd('set info on')
+ cubit.cmd('set echo on')
def surface_write(self,pathdir=None):
# optional surfaces, e.g. moho_surface
# should be created like e.g.:
@@ -756,9 +978,7 @@
faces=cubit.get_sub_elements('hex',h,2)
for f in faces:
if dic_quads_all.has_key(f):
- nodes=cubit.get_connectivity('face',f)
- txt='%10i %10i %10i %10i %10i\n' % (h,nodes[0],\
- nodes[1],nodes[2],nodes[3])
+ txt=self.create_facenode_string(h,f,cknormal=False)
surfhex_local.write(txt)
# closes file
surfhex_local.close()
@@ -781,7 +1001,10 @@
self.material_write(path+self.material_name)
self.nodescoord_write(path+self.nodecoord_name)
self.free_write(path+self.freename)
- self.abs_write(path+self.absname)
+ if self.cpml:
+ self.abs_write(path+self.cpmlname)
+ else:
+ self.abs_write(path+self.absname)
self.nummaterial_write(path+self.nummaterial_name)
# any other surfaces: ***surface***
self.surface_write(path)
@@ -789,17 +1012,17 @@
cubit.cmd('set info on')
cubit.cmd('set echo on')
-def export2SPECFEM3D(path_exporting_mesh_SPECFEM3D='.'):
- sem_mesh=mesh()
+def export2SPECFEM3D(path_exporting_mesh_SPECFEM3D='.',hex27=False,cpml=False,cpml_size=False,top_absorbing=False):
+ sem_mesh=mesh(hex27,cpml,cpml_size,top_absorbing)
#sem_mesh.block_definition()
#print sem_mesh.block_mat
#print sem_mesh.block_flag
#
sem_mesh.write(path=path_exporting_mesh_SPECFEM3D)
print 'END SPECFEM3D exporting process......'
-
+ if cpml:
+ cmd='save as "cpml.cub" overwrite'
+ cubit.cmd(cmd)
-if __name__ == '__main__':
- path='.'
- export2SPECFEM3D(path)
+
Added: seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/dev/cmpl.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/dev/cmpl.py (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/dev/cmpl.py 2013-09-26 00:57:45 UTC (rev 22850)
@@ -0,0 +1,198 @@
+#!/usr/bin/env python
+#############################################################################
+# cpml.py
+# this file is part of GEOCUBIT #
+# #
+# Created by Emanuele Casarotti #
+# Copyright (c) 2008 Istituto Nazionale di Geofisica e Vulcanologia #
+# #
+#############################################################################
+# #
+# GEOCUBIT is free software: you can redistribute it and/or modify #
+# it under the terms of the GNU General Public License as published by #
+# the Free Software Foundation, either version 3 of the License, or #
+# (at your option) any later version. #
+# #
+# GEOCUBIT is distributed in the hope that it will be useful, #
+# but WITHOUT ANY WARRANTY; without even the implied warranty of #
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
+# GNU General Public License for more details. #
+# #
+# You should have received a copy of the GNU General Public License #
+# along with GEOCUBIT. If not, see <http://www.gnu.org/licenses/>. #
+# #
+#############################################################################
+#
+try:
+ import start as start
+ cubit = start.start_cubit()
+except:
+ try:
+ import cubit
+ except:
+ print 'error importing cubit, check if cubit is installed'
+ pass
+
+
+def create_single_pml(n,coord,operator,limit,normal,distance,layers):
+ idv1=cubit.get_last_id('volume')+1
+ if normal=='0 0 -1':
+ txt="create element extrude face in surface with %s %s %s direction %s distance %s layers %s"
+ cmd=txt % (coord,operator,limit,normal,distance,layers)
+ else:
+ txt="create element extrude face in (surf in vol %s) with %s %s %s direction %s distance %s layers %s"
+ cmd=txt % (n,coord,operator,limit,normal,distance,layers)
+ cubit.cmd(cmd)
+ cmd='create mesh geometry Hex all except hex in vol all feature_angle 135.0'
+ cubit.cmd(cmd)
+ cmd='del hex all except hex in vol all'
+ cubit.cmd(cmd)
+ idv2=cubit.get_last_id('volume')
+ idv=' '.join(str(x) for x in range(idv1,idv2+1))
+ return idv
+
+def create_volume_pml(*args,**keys):
+ ymin=keys.get("ymin",False)
+ xmin=keys.get("xmin",False)
+ ymax=keys.get("ymax",False)
+ xmax=keys.get("xmax",False)
+ zmin=keys.get("zmin",False)
+ zmax=keys.get("zmax",False)
+ nvol=keys.get("vol",False)
+ layers=int(keys.get("layers",2))
+ size=float(keys.get("size",0))
+ thres=keys.get("thres",1)
+ distance=str(size*layers)
+ layers=str(layers)
+ #
+ #
+ for n in nvol:
+ normal='1 0 0'
+ coord='X_coord'
+ operator='>'
+ limit=str(xmax-thres)
+ idv=create_single_pml(n,coord,operator,limit,normal,distance,layers)
+ ##
+ ##
+ normal='-1 0 0'
+ coord='X_coord'
+ operator='<'
+ limit=str(xmin+thres)
+ idv2=create_single_pml(n,coord,operator,limit,normal,distance,layers)
+ #
+ #
+ #
+ normal='0 -1 0'
+ coord='Y_coord'
+ operator='<'
+ limit=str(ymin+thres)
+ idv=create_single_pml(n,coord,operator,limit,normal,distance,layers)
+ #
+ #
+ normal='0 1 0'
+ coord='Y_coord'
+ operator='>'
+ limit=str(ymax-thres)
+ idv=create_single_pml(n,coord,operator,limit,normal,distance,layers)
+ #
+ #
+ normal='0 0 -1'
+ coord='Z_coord'
+ operator='<'
+ limit=str(zmin+thres)
+ nvol='all'
+ idv=create_single_pml(nvol,coord,operator,limit,normal,distance,layers)
+
+
+def mesh_cpml(list_vol,remesh=True,refinement=None,top_surf=None,size=None):
+ from utilities import list2str
+ top_surf=list2str(top_surf)
+ if remesh:
+ cubit.cmd('reset vol all')
+ cubit.cmf('set dev on')
+ cubit.cmd('imprint vol all')
+ cubit.cmd('merge vol all')
+ cubit.cmd('vol all size '+str(size))
+ cubit.cmd('mesh vol all')
+ try:
+ for refdepth in refinement:
+ cubit.cmd('refine surf '+top_surf+' numsplit 1 bias 1 depth '+str(refdepth))
+ except:
+ print 'DEBUG: error in refinement cpml'
+ xmin=xmin-size
+ xmax=xmax+size
+ ymin=ymin-size
+ ymax=ymax+size
+ zmin=zmin-size
+ zmin=zmax-size
+ txt="group 'vol_xmin' add vol in surf with X_coord < "+str(xmin)
+ cubit.cmd(txt)
+ txt="group 'vol_xmax' add vol in surf with X_coord > "+str(xmax)
+ cubit.cmd(txt)
+ txt="group 'vol_ymin' add vol in surf with Y_coord < "+str(ymin)
+ cubit.cmd(txt)
+ txt="group 'vol_ymax' add vol in surf with Y_coord > "+str(ymax)
+ cubit.cmd(txt)
+ txt="group 'vol_zmin' add vol in surf with Z_coord < "+str(zmin)
+ cubit.cmd(txt)
+
+
+
+def collecting_cpml(ip,size=None,cpuxmin=0,cpuxmax=1,cpuymin=0,cpuymax=1,cpux=1,cpuy=1,cubfiles=False,decimate=False,layers=2):
+ import glob
+ import re
+ from utilities import list2str
+ #
+ if not size:
+ print 'cpml size must be specified'
+ return
+
+
+ boundary_dict={}
+ ##
+ try:
+ from boundary_definition import check_bc, map_boundary, define_surf
+ except:
+ pass
+ #
+ xmin,xmax,ymin,ymax,listfull=map_boundary(cpuxmin,cpuxmax,cpuymin,cpuymax,cpux,cpuy)
+ #
+ if cubfiles:
+ nf,listip,filenames,cubflag=importing_cubfiles(cubfiles)
+ else:
+ nf=0
+ filenames=[]
+ ip=0
+ #
+ if nf > 0:
+ for ip,filename in zip(listip,filenames):
+ try:
+ if ip in listfull:
+ if cubflag:
+ cubit.cmd('import cubit "'+filename+'"')
+ else:
+ cubit.cmd('import mesh geometry "'+filename+'" block all use nodeset sideset feature_angle 135.00 linear merge')
+ if decimate: cubit.cmd('refine volume all numsplit 1 bias 1.0 depth 1 ')
+ except:
+ cubit.cmd('import mesh geometry "'+filename+'" block all use nodeset sideset feature_angle 135.00 linear merge')
+ if decimate: cubit.cmd('refine volume all numsplit 1 bias 1.0 depth 1 ')
+ ip=0
+ if decimate: cubit.cmd('export mesh "decimated_before_cmpl.e" dimension 3 block all overwrite')
+ else:
+ if decimate:
+ cubit.cmd('refine volume all numsplit 1 bias 1.0 depth 1 ')
+ cubit.cmd('export mesh "decimated_before_cmpl.e" dimension 3 block all overwrite')
+
+
+ #
+ #
+ #print boundary_dict
+ block_list=cubit.get_block_id_list()
+ for block in block_list:
+ ty=cubit.get_block_element_type(block)
+ if ty == 'HEX8':
+ cubit.cmd('block '+str(block)+' name "vol'+str(block)+'"')
+
+ list_vol=list(cubit.parse_cubit_list('volume','all'))
+ create_pml(xmin=xmin,xmax=xmax,ymax=ymax,ymin=ymin,zmin=zmin,zmax=zmax,size=size,layers=layers,vol=list_vol)
+
Property changes on: seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/dev/cmpl.py
___________________________________________________________________
Added: svn:executable
+ *
Deleted: seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/dev/cubit2googleearth.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/dev/cubit2googleearth.py 2013-09-26 00:49:25 UTC (rev 22849)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/dev/cubit2googleearth.py 2013-09-26 00:57:45 UTC (rev 22850)
@@ -1,297 +0,0 @@
-#!python
-#############################################################################
-# cubit2specfem3d.py #
-# this file is part of GEOCUBIT #
-# #
-# Created by Emanuele Casarotti #
-# Copyright (c) 2008 Istituto Nazionale di Geofisica e Vulcanologia #
-# #
-#############################################################################
-# #
-# GEOCUBIT is free software: you can redistribute it and/or modify #
-# it under the terms of the GNU General Public License as published by #
-# the Free Software Foundation, either version 3 of the License, or #
-# (at your option) any later version. #
-# #
-# GEOCUBIT is distributed in the hope that it will be useful, #
-# but WITHOUT ANY WARRANTY; without even the implied warranty of #
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
-# GNU General Public License for more details. #
-# #
-# You should have received a copy of the GNU General Public License #
-# along with GEOCUBIT. If not, see <http://www.gnu.org/licenses/>. #
-# #
-#############################################################################
-#
-#
-try:
- cubit.cmd('comment')
-except:
- import cubit
-
-KMLstart='<kml xmlns="http://earth.google.com/kml/2.2">\n<Folder><name>prova</name>\n'
-placemarkstart='<Placemark>\n <name></name>\n <description></description>\n<TimeStamp><begin>%s</begin><end>%s</end></TimeStamp>\n <MultiGeometry>\n <altitudeMode>absolute</altitudeMode>\n'
-placemarkend=' </MultiGeometry>\n</Placemark>\n'
-templatepoly=' <Polygon>\n <altitudeMode>absolute</altitudeMode>\n <outerBoundaryIs>\n <LinearRing>\n <coordinates>\n %s,%s,%s\n %s,%s,%s\n %s,%s,%s\n %s,%s,%s\n %s,%s,%s \n </coordinates>\n </LinearRing>\n </outerBoundaryIs>\n </Polygon>\n'
-
-KMLend=' </Folder>\n</kml>\n'
-
-
-import datetime
-now=datetime.datetime.now().year
-
-import random
-
-# Lat Long - UTM, UTM - Lat Long conversions
-
-from math import pi, sin, cos, tan, sqrt
-
-#LatLong- UTM conversion..h
-#definitions for lat/long to UTM and UTM to lat/lng conversions
-#include <string.h>
-
-_deg2rad = pi / 180.0
-_rad2deg = 180.0 / pi
-
-_EquatorialRadius = 2
-_eccentricitySquared = 3
-
-_ellipsoid = [[ -1, "Placeholder", 0, 0],[ 1, "Airy", 6377563, 0.00667054],[ 2, "Australian National", 6378160, 0.006694542],[ 3, "Bessel 1841", 6377397, 0.006674372],[ 4, "Bessel 1841 (Nambia] ", 6377484, 0.006674372],[ 5, "Clarke 1866", 6378206, 0.006768658],[ 6, "Clarke 1880", 6378249, 0.006803511],[ 7, "Everest", 6377276, 0.006637847],[ 8, "Fischer 1960 (Mercury] ", 6378166, 0.006693422],[ 9, "Fischer 1968", 6378150, 0.006693422],[ 10, "GRS 1967", 6378160, 0.006694605],[ 11, "GRS 1980", 6378137, 0.00669438],[ 12, "Helmert 1906", 6378200, 0.006693422],[ 13, "Hough", 6378270, 0.00672267],[ 14, "International", 6378388, 0.00672267],[ 15, "Krassovsky", 6378245, 0.006693422],[ 16, "Modified Airy", 6377340, 0.00667054],[ 17, "Modified Everest", 6377304, 0.006637847],[ 18, "Modified Fischer 1960", 6378155, 0.006693422],[ 19, "South American 1969", 6378160, 0.006694542],[ 20, "WGS 60", 6378165, 0.006693422],[ 21, "WGS 66", 6378145, 0.006694542],[ 22, "WGS-72", 6378135, 0.006694318],[ 23, "WGS-84", 6378137, 0.00669438]]
-
-#Reference ellipsoids derived from Peter H. Dana's website-
-#http://www.utexas.edu/depts/grg/gcraft/notes/datum/elist.htmlV
-#Department of Geography, University of Texas at Austin
-#Internet: pdana at mail.utexas.edu
-#3/22/95
-
-#Source
-#Defense Mapping Agency. 1987b. DMA Technical Report: Supplement to Department of Defense World Geodetic System
-#1984 Technical Report. Part I and II. Washington, DC: Defense Mapping Agency
-
-def LLtoUTM(ReferenceEllipsoid, Lat, Long):
- #converts lat/long to UTM coords. Equations from USGS Bulletin 1532
- #East Longitudes are positive, West longitudes are negative.
- #North latitudes are positive, South latitudes are negative
- #Lat and Long are in decimal degrees
- #Written by Chuck Gantz- chuck.gantz at globalstar.com
- a = _ellipsoid[ReferenceEllipsoid][_EquatorialRadius]
- eccSquared = _ellipsoid[ReferenceEllipsoid][_eccentricitySquared]
- k0 = 0.9996#Make sure the longitude is between -180.00 .. 179.9
- LongTemp = (Long+180)-int((Long+180)/360)*360-180 # -180.00 .. 179.9
- LatRad = Lat*_deg2rad
- LongRad = LongTemp*_deg2rad
- ZoneNumber = int((LongTemp + 180)/6) + 1
- if Lat >= 56.0 and Lat < 64.0 and LongTemp >= 3.0 and LongTemp < 12.0:
- ZoneNumber = 32
- # Special zones for Svalbard
- if Lat >= 72.0 and Lat < 84.0:
- if LongTemp >= 0.0 and LongTemp < 9.0:ZoneNumber = 31
- elif LongTemp >= 9.0 and LongTemp < 21.0: ZoneNumber = 33
- elif LongTemp >= 21.0 and LongTemp < 33.0: ZoneNumber = 35
- elif LongTemp >= 33.0 and LongTemp < 42.0: ZoneNumber = 37
- LongOrigin = (ZoneNumber - 1)*6 - 180 + 3 #+3 puts origin in middle of zone
- LongOriginRad = LongOrigin * _deg2rad
- #compute the UTM Zone from the latitude and longitude
- UTMZone = "%d%c" % (ZoneNumber, _UTMLetterDesignator(Lat))
- eccPrimeSquared = (eccSquared)/(1-eccSquared)
- N = a/sqrt(1-eccSquared*sin(LatRad)*sin(LatRad))
- T = tan(LatRad)*tan(LatRad)
- C = eccPrimeSquared*cos(LatRad)*cos(LatRad)
- A = cos(LatRad)*(LongRad-LongOriginRad)
- M = a*((1
- - eccSquared/4
- - 3*eccSquared*eccSquared/64
- - 5*eccSquared*eccSquared*eccSquared/256)*LatRad
- - (3*eccSquared/8
- + 3*eccSquared*eccSquared/32
- + 45*eccSquared*eccSquared*eccSquared/1024)*sin(2*LatRad)
- + (15*eccSquared*eccSquared/256 + 45*eccSquared*eccSquared*eccSquared/1024)*sin(4*LatRad)
- - (35*eccSquared*eccSquared*eccSquared/3072)*sin(6*LatRad))
- UTMEasting = (k0*N*(A+(1-T+C)*A*A*A/6
- + (5-18*T+T*T+72*C-58*eccPrimeSquared)*A*A*A*A*A/120)
- + 500000.0)
- UTMNorthing = (k0*(M+N*tan(LatRad)*(A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24
- + (61
- -58*T
- +T*T
- +600*C
- -330*eccPrimeSquared)*A*A*A*A*A*A/720)))
- if Lat < 0:
- UTMNorthing = UTMNorthing + 10000000.0; #10000000 meter offset for southern hemisphere
- return (UTMZone, UTMEasting, UTMNorthing)
-
-
-def _UTMLetterDesignator(Lat):
-#This routine determines the correct UTM letter designator for the given latitude
-#returns 'Z' if latitude is outside the UTM limits of 84N to 80S
-#Written by Chuck Gantz- chuck.gantz at globalstar.com
- if 84 >= Lat >= 72: return 'X'
- elif 72 > Lat >= 64: return 'W'
- elif 64 > Lat >= 56: return 'V'
- elif 56 > Lat >= 48: return 'U'
- elif 48 > Lat >= 40: return 'T'
- elif 40 > Lat >= 32: return 'S'
- elif 32 > Lat >= 24: return 'R'
- elif 24 > Lat >= 16: return 'Q'
- elif 16 > Lat >= 8: return 'P'
- elif 8 > Lat >= 0: return 'N'
- elif 0 > Lat >= -8: return 'M'
- elif -8> Lat >= -16: return 'L'
- elif -16 > Lat >= -24: return 'K'
- elif -24 > Lat >= -32: return 'J'
- elif -32 > Lat >= -40: return 'H'
- elif -40 > Lat >= -48: return 'G'
- elif -48 > Lat >= -56: return 'F'
- elif -56 > Lat >= -64: return 'E'
- elif -64 > Lat >= -72: return 'D'
- elif -72 > Lat >= -80: return 'C'
- else: return 'Z' # if the Latitude is outside the UTM limits
-
-#void UTMtoLL(int ReferenceEllipsoid, const double UTMNorthing, const double UTMEasting, const char* UTMZone,
-# double& Lat, double& Long )
-
-def UTMtoLL(ReferenceEllipsoid, northing, easting, zone):
- #converts UTM coords to lat/long. Equations from USGS Bulletin 1532
- #East Longitudes are positive, West longitudes are negative.
- #North latitudes are positive, South latitudes are negative
- #Lat and Long are in decimal degrees.
- #Written by Chuck Gantz- chuck.gantz at globalstar.com
- #Converted to Python by Russ Nelson <nelson at crynwr.com>
- k0 = 0.9996
- a = _ellipsoid[ReferenceEllipsoid][_EquatorialRadius]
- eccSquared = _ellipsoid[ReferenceEllipsoid][_eccentricitySquared]
- e1 = (1-sqrt(1-eccSquared))/(1+sqrt(1-eccSquared))
- #NorthernHemisphere; //1 for northern hemispher, 0 for southern
- x = easting - 500000.0 #remove 500,000 meter offset for longitude
- y = northing
- ZoneLetter = zone[-1]
- ZoneNumber = int(zone[:-1])
- if ZoneLetter >= 'N':
- NorthernHemisphere = 1 # point is in northern hemisphere
- else:
- NorthernHemisphere = 0 # point is in southern hemisphere
- y -= 10000000.0 # remove 10,000,000 meter offset used for southern hemisphere
- LongOrigin = (ZoneNumber - 1)*6 - 180 + 3 # +3 puts origin in middle of zone
- eccPrimeSquared = (eccSquared)/(1-eccSquared)
- M = y / k0
- mu = M/(a*(1-eccSquared/4-3*eccSquared*eccSquared/64-5*eccSquared*eccSquared*eccSquared/256))
- phi1Rad = (mu + (3*e1/2-27*e1*e1*e1/32)*sin(2*mu)
- + (21*e1*e1/16-55*e1*e1*e1*e1/32)*sin(4*mu)
- +(151*e1*e1*e1/96)*sin(6*mu))
- phi1 = phi1Rad*_rad2deg;
- N1 = a/sqrt(1-eccSquared*sin(phi1Rad)*sin(phi1Rad))
- T1 = tan(phi1Rad)*tan(phi1Rad)
- C1 = eccPrimeSquared*cos(phi1Rad)*cos(phi1Rad)
- R1 = a*(1-eccSquared)/pow(1-eccSquared*sin(phi1Rad)*sin(phi1Rad), 1.5)
- D = x/(N1*k0)
- Lat = phi1Rad - (N1*tan(phi1Rad)/R1)*(D*D/2-(5+3*T1+10*C1-4*C1*C1-9*eccPrimeSquared)*D*D*D*D/24
- +(61+90*T1+298*C1+45*T1*T1-252*eccPrimeSquared-3*C1*C1)*D*D*D*D*D*D/720)
- Lat = Lat * _rad2deg
- Long = (D-(1+2*T1+C1)*D*D*D/6+(5-2*C1+28*T1-3*C1*C1+8*eccPrimeSquared+24*T1*T1)
- *D*D*D*D*D/120)/cos(phi1Rad)
- Long = LongOrigin + Long * _rad2deg
- return (Lat, Long)
-
-class cubit2GE(object):
- def __init__(self,unit='utm',x_origin=None,y_origin=None,z_origin=None):
- super(cubit2GE, self).__init__()
- self.unit=unit
- self.xo=x_origin
- self.yo=y_origin
- self.zo=z_origin
- cubit.cmd('compress all')
- def __repr__(self):
- pass
- def create_blocks(self,list_entity=None):
- txt=' face in surface '
- if list_entity:
- if not isinstance(list_entity,list):
- list_entity=[list_entity]
- iblock=cubit.get_next_block_id()
- for entity in list_entity:
- command = "block "+str(iblock)+ txt +str(entity)
- cubit.cmd(command)
- command = "block "+str(iblock)+ "name 'toGE'"
- cubit.cmd(command)
- self.iblock=iblock
- def geo2utm(self,lon,lat,unit,ellipsoid=23):
- if unit == 'geo' :
- (zone, x, y) = LLtoUTM(ellipsoid, lat, lon)
- elif unit == 'utm' :
- x=lon
- y=lat
- zone=0
- return zone,x,y
- def utm2geo(self,lon,lat,zone,unit='utm'):
- if unit == 'utm' and zone:
- (y, x) = UTMtoLL(23, lat, lon,zone)
- else:
- x=lon
- y=lat
- return x,y
- def write(self,GEname=None,t1=0,t2=1):
- cubit.cmd('set info off')
- cubit.cmd('set echo off')
- cubit.cmd('set journal off')
- from sets import Set
- if not GEname: GEname='cubit2GE.kml'
- if self.xo and self.yo:
- zone_utm,xo_utm,yo_utm=self.geo2utm(self.xo,self.yo,'geo')
- else:
- xo_utm=0.
- yo_utm=0.
- zone_utm=None
- if self.zo:
- zo=self.zo
- else:
- zo=0
- #
- GE=open(GEname,'w')
- print 'Writing '+GEname+'.....'
- #
- #
- GE.write(KMLstart)
- for it in range(int(t1),int(t2),1):
- tbegin=now+seconds*it
- tafter=now+(seconds*(it+1))
- print tbegin.isoformat()
- GE.write(placemarkstart%(str(tbegin.isoformat()),str(tafter.isoformat())))
- #
- block=self.iblock
- nameb=cubit.get_exodus_entity_name('block',block)
- print nameb,block
- quads_all=cubit.get_block_faces(block)
- print len(quads_all)
- for f in quads_all:
- nodes=cubit.get_connectivity('Face',f)
- coord=[]
- for n in [nodes[0],nodes[1],nodes[2],nodes[3],nodes[0]]:
- x,y,z=cubit.get_nodal_coordinates(n)
- if self.unit == 'geo':
- zone,x,y=self.geo2utm(x,y,'geo')
- new_x=x+xo_utm+sin(it)*3
- new_y=y+yo_utm+cos(it)*3
- new_z=z+zo+sin(it)*3
- lon,lat=self.utm2geo(new_x,new_y,zone_utm,unit='utm')
- coord.append(str(lon))
- coord.append(str(lat))
- coord.append(str(new_z))
- GE.write(templatepoly % tuple(coord))
- GE.write(placemarkend)
- GE.write(KMLend)
- GE.close()
- #
- print 'Ok'
- cubit.cmd('set info on')
- cubit.cmd('set echo on')
-
-def toGE(name='c2GE.kml',list_entity=None,unit='utm',x_origin=None,y_origin=None,z_origin=None,t1=0,t2=5):
- GE=cubit2GE(unit=unit,x_origin=x_origin,y_origin=y_origin,z_origin=z_origin)
- GE.create_blocks(list_entity=list_entity)
- GE.write(GEname=name,t1=t1,t2=t2)
-
-name='/Users/emanuele/Desktop/c2GE.kml'
-toGE(name,list_entity='1 to 493',x_origin=8.767022,y_origin= 44.460383, z_origin=220)
-
-
-
Deleted: seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/dev/material.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/dev/material.py 2013-09-26 00:49:25 UTC (rev 22849)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/dev/material.py 2013-09-26 00:57:45 UTC (rev 22850)
@@ -1,80 +0,0 @@
-#############################################################################
-# material.py
-# this file is part of GEOCUBIT #
-# #
-# Created by Emanuele Casarotti #
-# Copyright (c) 2008 Istituto Nazionale di Geofisica e Vulcanologia #
-# #
-#############################################################################
-# #
-# GEOCUBIT is free software: you can redistribute it and/or modify #
-# it under the terms of the GNU General Public License as published by #
-# the Free Software Foundation, either version 3 of the License, or #
-# (at your option) any later version. #
-# #
-# GEOCUBIT is distributed in the hope that it will be useful, #
-# but WITHOUT ANY WARRANTY; without even the implied warranty of #
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
-# GNU General Public License for more details. #
-# #
-# You should have received a copy of the GNU General Public License #
-# along with GEOCUBIT. If not, see <http://www.gnu.org/licenses/>. #
-# #
-#############################################################################
-
-
-
-class material():
- def __init__(self,name,vp,vs,rho,q,ani):
- self.name=name
- self.vp=vp
- self.vs=vs
- self.rho=rho
- self.qp=q
- self.qs=ani
- def __repr__(self):
- print self.name
- print 'vp',self.vp
- print 'vs',self.vs
- print 'rho',self.rho
- print 'Qp',self.qp
- print 'Qs',self.qs
- return
- def __str__(self):
- txt=self.name+' vp '+str(self.vp)+ ' vs '+str(self.vs)+ ' rho '+str(self.rho)+ ' Qp '+str(self.qp)+ ' Qs '+str(self.qs)
- return txt
-
-
-def read_material_archive(archivefile='material_archive.dat'):
- archive=open(archivefile,'r')
- archive_dict={}
- for record in archive:
- m=material(*record.split())
- archive_dict[m.name]=m
- return archive_dict
-
-
-def block_description(archivefile='material_archive.dat'):
- try:
- import initializing as initializing
- cubit = initializing.initializing_cubit()
- except:
- pass
- import menu as menu
- archive=read_material_archive(archivefile=menu.material_file)
- b_assign=menu.material_assignement
- print b_assign
- if len(b_assign) != 0:
- for block in b_assign:
- mat=archive[block[1]]
- print mat
- cubit.cmd("block "+block[0]+" attribute count 6")
- cubit.cmd("block "+block[0]+" name '"+ mat.name +'"')
- cubit.cmd("block "+block[0]+" attribute index 1 "+ '1' )
- cubit.cmd("block "+block[0]+" attribute index 2 "+ str(mat.vp) )
- cubit.cmd("block "+block[0]+" attribute index 3 "+ str(mat.vs) )
- cubit.cmd("block "+block[0]+" attribute index 4 "+ str(mat.rho) )
- cubit.cmd("block "+block[0]+" attribute index 5 "+ str(mat.qp) )
- cubit.cmd("block "+block[0]+" attribute index 5 "+ str(mat.qs) )
-
-
\ No newline at end of file
Modified: seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/exportlib.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/exportlib.py 2013-09-26 00:49:25 UTC (rev 22849)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/exportlib.py 2013-09-26 00:57:45 UTC (rev 22850)
@@ -33,6 +33,86 @@
print 'error importing cubit, check if cubit is installed'
pass
+import glob
+
+def add_sea_layer(block=1001,optionsea=False):
+ if optionsea:
+ sea=optionsea['sea']
+ seaup=optionsea['seaup']
+ sealevel=optionsea['sealevel']
+ seathres=optionsea['seathres']
+ else:
+ sea=False
+ seaup=False
+ sealevel=False
+ seathres=False
+
+ ######TODO
+ #add sea hex
+ #change hex absoorbing....
+ block_list=cubit.get_block_id_list()
+ id_block = max(block for block in block_list if block<1000)
+ cubit.cmd('delete block '+str(id_block))
+ #sea
+ command= 'block '+str(id_block)+' hex in node in face in block '+str(block)+' with Z_coord < '+str(seathres)
+ cubit.cmd(command)
+ command = "block "+str(id_block)+" name 'sea'"
+ cubit.cmd(command)
+ if not seaup:
+ id_block+=1
+ command= 'block '+str(id_block)+' hex in node in face in block '+str(block)+' with (Z_coord > '+str(seathres)+' and Z_coord < '+str(sealevel)+')'
+ cubit.cmd(command)
+ command = "block "+str(id_block)+" name 'shwater'"
+ cubit.cmd(command)
+ id_block+=1
+ command= 'block '+str(id_block)+' hex in node in face in block '+str(block)+' with Z_coord >= '+str(sealevel)
+ cubit.cmd(command)
+ command = "block "+str(id_block)+" name 'continent'"
+ cubit.cmd(command)
+
+
+def importing_cubfiles(cubfiles):
+ import re
+ rule_st=re.compile("(.+)_[0-9]+\.")
+ rule_ex=re.compile(".+_[0-9]+\.(.+)")
+ rule_int=re.compile(".+_([0-9]+)\.")
+ filenames=glob.glob(cubfiles)
+ try:
+ st = rule_st.findall(filenames[0])[0]
+ ex = rule_ex.findall(filenames[0])[0]
+ listflag=True
+ except:
+ ex=''
+ listflag=False
+ if ex == 'cub':
+ cubflag=True
+ else:
+ cubflag=False
+ list_int=[]
+ fs=[]
+ try:
+ for f in filenames:
+ i=int(rule_int.findall(f)[0])
+ list_int.append(i)
+ list_int.sort()
+ for i,ind in enumerate(list_int):
+ f=st+'_'+str(ind)+'.'+ex
+ fs.append(f)
+ except:
+ pass
+ if listflag:
+ filenames=fs
+ else:
+ pass
+ return len(filenames),list_int,filenames,cubflag
+
+
+
+
+
+
+
+
def refine_closecurve(block=1001,closed_filenames=None,acis=True):
from utilities import load_curves
from boundary_definition import build_block_side,define_surf
@@ -150,7 +230,7 @@
-def collecting_merging(cpuxmin=0,cpuxmax=1,cpuymin=0,cpuymax=1,cpux=1,cpuy=1,cubfiles=False,ckbound_method1=False,ckbound_method2=False,merge_tolerance=None):
+def collecting_merging(cpuxmin=0,cpuxmax=1,cpuymin=0,cpuymax=1,cpux=1,cpuy=1,cubfiles=False,ckbound_method1=False,ckbound_method2=False,merge_tolerance=None,decimate=False):
import glob
import re
#
@@ -167,49 +247,21 @@
xmin,xmax,ymin,ymax,listfull=map_boundary(cpuxmin,cpuxmax,cpuymin,cpuymax,cpux,cpuy)
#
if cubfiles:
- filenames=glob.glob(cubfiles)
- try:
- st = rule_st.findall(filenames[0])[0]
- ex = rule_ex.findall(filenames[0])[0]
- listflag=True
- except:
- ex=''
- listflag=False
-
- if ex == 'cub':
- cubflag=True
- else:
- cubflag=False
- list_int=[]
- fs=[]
- try:
- for f in filenames:
- i=int(rule_int.findall(f)[0])
- list_int.append(i)
- list_int.sort()
- for i,ind in enumerate(list_int):
- f=st+'_'+str(ind)+'.'+ex
- fs.append(f)
- except:
- pass
- if listflag:
- filenames=fs
- else:
- pass
- print len(filenames),filenames
+ nf,listip,filenames,cubflag=importing_cubfiles(cubfiles)
else:
+ nf=0
filenames=[]
ip=0
#
- if len(filenames) > 0:
- for filename in filenames[:]:
+ if nf > 0:
+ for ip,filename in zip(listip,filenames):
try:
- ip=int(rule_int.findall(filename)[0])
if ip in listfull:
if cubflag:
cubit.cmd('import cubit "'+filename+'"')
else:
cubit.cmd('import mesh geometry "'+filename+'" block all use nodeset sideset feature_angle 135.00 linear merge')
+ if decimate: cubit.cmd('refine volume all numsplit 1 bias 1.0 depth 1 ')
boundary=check_bc(ip,xmin,xmax,ymin,ymax,cpux,cpuy,cpuxmin,cpuxmax,cpuymin,cpuymax)
boundary_dict[ip]=boundary
list_vol=list(cubit.parse_cubit_list('volume','all'))
@@ -219,6 +271,7 @@
cubit.cmd(command)
except:
cubit.cmd('import mesh geometry "'+filename+'" block all use nodeset sideset feature_angle 135.00 linear merge')
+ if decimate: cubit.cmd('refine volume all numsplit 1 bias 1.0 depth 1 ')
ip=0
boundary=check_bc(ip,xmin,xmax,ymin,ymax,cpux,cpuy,cpuxmin,cpuxmax,cpuymin,cpuymax)
boundary_dict[ip]=boundary
@@ -229,9 +282,11 @@
cubit.cmd(command)
cubit.cmd('export mesh "tmp_collect_NOmerging.e" dimension 3 block all overwrite')
else:
+ if decimate: cubit.cmd('refine volume all numsplit 1 bias 1.0 depth 1 ')
boundary=check_bc(ip,xmin,xmax,ymin,ymax,cpux,cpuy,cpuxmin,cpuxmax,cpuymin,cpuymax)
#
#
+ #print boundary_dict
block_list=cubit.get_block_id_list()
for block in block_list:
ty=cubit.get_block_element_type(block)
@@ -240,7 +295,10 @@
#
#
print 'chbound',ckbound_method1,ckbound_method2
- if ckbound_method1 and not ckbound_method2:
+
+
+ if ckbound_method1 and not ckbound_method2 and len(filenames) != 1:
+ #use the equivalence method for groups
if isinstance(merge_tolerance,list):
tol=merge_tolerance[0]
elif merge_tolerance:
@@ -281,11 +339,12 @@
idown=ip
idiag=ileft
#
+ print ip,ileft,idiag,idown
if ip != idown:
nup=boundary_dict[ip]['nodes_surf_ymin']
ndow=boundary_dict[idown]['nodes_surf_ymax']
- merge_node(nup,ndow)
-
+ merge_node_ck(nup,ndow)
+
if idiag != idown:
if ip in ymax and ip not in xmin:
nlu=boundary_dict[ip]['node_curve_xminymax'] #node in curve chunck left up... r u
@@ -301,15 +360,21 @@
nlu=boundary_dict[ileft]['node_curve_xmaxymin']
merge_node_4(nru,nrd,nld,nlu)
elif ip in xmin:
- nru=boundary_dict[ip]['node_curve_xminymin'] #node in curve chunck right up... r u
- nrd=boundary_dict[idown]['node_curve_xminymax']
+ nlu=boundary_dict[ip]['node_curve_xminymin'] #node in curve chunck right up... r u
+ nld=boundary_dict[idown]['node_curve_xminymax']
+ merge_node(nld,nlu)
+ nru=boundary_dict[ip]['node_curve_xmaxymin'] #node in curve chunck right up... r u
+ nrd=boundary_dict[idown]['node_curve_xmaxymax']
merge_node(nrd,nru)
+
+
+
#
if ip != ileft:
nright=boundary_dict[ip]['nodes_surf_xmin']
nleft=boundary_dict[ileft]['nodes_surf_xmax']
- merge_node(nright,nleft)
+ merge_node_ck(nright,nleft)
#
#
if ip in ymin:
@@ -334,7 +399,7 @@
n1=cubit.get_group_nodes(group_id_1)
if len(n1) != 0:
print 'error, negative jacobian after the equivalence node command, use --merge2 instead of --equivalence/--merge/--merge1'
- elif ckbound_method2 and not ckbound_method1:
+ elif ckbound_method2 and not ckbound_method1 and len(filenames) != 1:
if isinstance(merge_tolerance,list):
tol=merge_tolerance[0]
elif merge_tolerance:
@@ -428,7 +493,7 @@
n1=cubit.get_group_nodes(group_id_1)
if len(n1) != 0:
print 'error, negative jacobian after the equivalence node command, check the mesh'
- elif ckbound_method1 and ckbound_method2:
+ elif ckbound_method1 and ckbound_method2 and len(filenames) != 1:
block_list=cubit.get_block_id_list()
i=-1
for block in block_list:
@@ -474,15 +539,27 @@
if len(n1) != 0:
print 'error, negative jacobian after the equivalence node command, use --merge instead of --equivalence'
-def collect(cpuxmin=0,cpuxmax=1,cpuymin=0,cpuymax=1,cpux=1,cpuy=1,cubfiles=False,ckbound_method1=False,ckbound_method2=False,merge_tolerance=None,curverefining=False,outfilename='totalmesh_merged',qlog=False,export2SPECFEM3D=False,listblock=None,listflag=None,outdir='.'):
+
+
+
+def collect(cpuxmin=0,cpuxmax=1,cpuymin=0,cpuymax=1,cpux=1,cpuy=1,cubfiles=False,ckbound_method1=False,ckbound_method2=False,merge_tolerance=None,curverefining=False,outfilename='totalmesh_merged',qlog=False,export2SPECFEM3D=False,listblock=None,listflag=None,outdir='.',add_sea=False,decimate=False,cpml=False,cpml_size=False,top_absorbing=False,hex27=False):
#
- collecting_merging(cpuxmin,cpuxmax,cpuymin,cpuymax,cpux,cpuy,cubfiles=cubfiles,ckbound_method1=ckbound_method1,ckbound_method2=ckbound_method2,merge_tolerance=merge_tolerance)
+ cubit.cmd('set journal error off')
+ cubit.cmd('set verbose error off')
+ collecting_merging(cpuxmin,cpuxmax,cpuymin,cpuymax,cpux,cpuy,cubfiles=cubfiles,ckbound_method1=ckbound_method1,ckbound_method2=ckbound_method2,merge_tolerance=merge_tolerance,decimate=decimate)
+ cubit.cmd('set journal error on')
+ cubit.cmd('set verbose error on')
#
if curverefining:
block=1001 #topography
refine_closecurve(block,curverefining,acis=True)
#
#
+ if add_sea:
+ block=1001
+ add_sea_layer(block=block)
+
+
cubit.cmd('compress all')
command="export mesh '"+outfilename+".e' block all overwrite xml '"+outfilename+".xml'"
cubit.cmd(command)
@@ -521,9 +598,9 @@
#
#
if export2SPECFEM3D:
- e2SEM(files=False,listblock=listblock,listflag=listflag,outdir=outdir)
+ e2SEM(files=False,listblock=listblock,listflag=listflag,outdir=outdir,cpml=cpml,cpml_size=cpml_size,top_absorbing=top_absorbing,hex27=hex27)
-def e2SEM(files=False,listblock=None,listflag=None,outdir='.'):
+def e2SEM(files=False,listblock=None,listflag=None,outdir='.',cpml=False,cpml_size=False,top_absorbing=False,hex27=False):
import glob
if files:
filenames=glob.glob(files)
@@ -547,14 +624,14 @@
if 'HEX' in ty:
listblock.append(block)
#listflag.append(block)
- listflag=range(1,len(listflag)+1)
+ listflag=range(1,len(block_list)+1)
#
for ib,iflag in zip(listblock,listflag):
cubit.cmd("block "+str(ib)+" attribute count 1")
cubit.cmd("block "+str(ib)+" attribute index 1 "+ str(iflag) )
#
import cubit2specfem3d
- cubit2specfem3d.export2SPECFEM3D(outdir)
+ cubit2specfem3d.export2SPECFEM3D(outdir,cpml=cpml,cpml_size=cpml_size,top_absorbing=top_absorbing,hex27=hex27)
def invert_dict(d):
inv = {}
@@ -600,38 +677,106 @@
return factor,minvalue,inv_length
-def merge_node(n1,n2):
+def merge_node_ck(n1,n2):
factor,minvalue,inv_length=prepare_equivalence(n1,n2)
+ cubit.cmd('set info off')
+ cubit.cmd('set echo off')
+ cubit.cmd('set journal off')
+ cubit.cmd('set error off')
+
for k in inv_length.keys()[:-1]:
if len(inv_length[k]) > 0:
cmd='equivalence node '+' '.join(' '.join(str(n) for n in x) for x in inv_length[k])+' tolerance '+str(k*factor+minvalue/2.)
cubit.cmd(cmd)
print 'equivalence '+str(len(inv_length[k]))+' couples of nodes - tolerance '+str(k*factor+minvalue/2.)
+
+ cubit.cmd('group "checkmerge" add node '+' '.join(str(n) for n in n1)+' '+' '.join(str(n) for n in n2))
+ idg=cubit.get_id_from_name('checkmerge')
+ remainnodes=cubit.get_group_nodes(idg)
+ print 'from '+str(len(n1)+len(n2))+' nodes -> '+str(len(remainnodes)) +' nodes'
+ if len(n1) != len(remainnodes):
+ print 'equivalence '+str(len(remainnodes))+' couples of nodes - tolerance '+str(minvalue/2.)
+ cubit.cmd('set info on')
+ cubit.cmd('set echo on')
+ cubit.cmd('set journal on')
+ cmd='equivalence node in group '+str(idg)+' tolerance '+str(minvalue/2.)
+ cubit.cmd(cmd)
+ cmd='block 3000 node in group '+str(idg)
+ cubit.cmd(cmd)
+
+ if len(n1) != len(remainnodes):
+ cubit.cmd('export mesh "error_merging.e" dimension 3 block all overwrite')
+ cubit.cmd('save as "error_merging.cub" dimension 3 block all overwrite')
+ print 'error merging '
+ if False:
+ import sys
+ sys.exit(2)
+
+ cubit.cmd('delete group checkmerge')
+ cubit.cmd('delete block 3000')
+
+ cubit.cmd('set info on')
+ cubit.cmd('set echo on')
+ cubit.cmd('set journal on')
+
+
+
+def merge_node(n1,n2):
+ factor,minvalue,inv_length=prepare_equivalence(n1,n2)
+
+ cubit.cmd('set info off')
+ cubit.cmd('set echo off')
+ cubit.cmd('set journal off')
+
+
+ for k in inv_length.keys()[:-1]:
+ if len(inv_length[k]) > 0:
+ cmd='equivalence node '+' '.join(' '.join(str(n) for n in x) for x in inv_length[k])+' tolerance '+str(k*factor+minvalue/2.)
+ cubit.cmd(cmd)
+ print 'equivalence '+str(len(inv_length[k]))+' couples of nodes - tolerance '+str(k*factor+minvalue/2.)
+
+ cubit.cmd('set info on')
+ cubit.cmd('set echo on')
+ cubit.cmd('set journal on')
+
+
+
+
+
def prepare_equivalence_4(nodes1,nodes2,nodes3,nodes4):
cubit.cmd('set info off')
cubit.cmd('set echo off')
cubit.cmd('set journal off')
length={}
- for ns in zip(nodes1,nodes2,nodes3,nodes4):
- cmd='group "tmpn" add edge in node '+' '.join(str(n) for n in ns )
- cubit.cmd(cmd)
- ge=cubit.get_id_from_name("tmpn")
- e1=cubit.get_group_edges(ge)
- lengthmin=1e9
- for e in e1:
- lengthmin=min(lengthmin,cubit.get_mesh_edge_length(e))
- length[ns]=lengthmin*.5
- cubit.cmd('delete group '+str(ge))
+ nodes=[nodes1,nodes2,nodes3,nodes4]
+ check=map(len,nodes)
+ checked_nodes=[]
+ for ind,iflag in enumerate(check):
+ if iflag:
+ checked_nodes=checked_nodes+nodes[ind]
+
+ cmd='group "tmpn" add edge in node '+' '.join(str(n) for n in checked_nodes )
+ cubit.cmd(cmd)
+ ge=cubit.get_id_from_name("tmpn")
+ e1=cubit.get_group_edges(ge)
+ lengthmin=1e9
+ for e in e1:
+ lengthmin=min(lengthmin,cubit.get_mesh_edge_length(e))
+ length[e]=lengthmin*.5
+ cubit.cmd('delete group '+str(ge))
try:
minvalue=min(length.values())
maxvalue=max(length.values())
except:
- print nodes1,nodes2,nodes3,nodes4
- print 'edges ', e1
- minvalue=100.
+ try:
+ print nodes
+ print 'edges ', e1
+ except:
+ pass
+ minvalue=10.
maxvalue=2000.
print 'min lentgh: ',minvalue,'max lentgh: ',maxvalue
nbin= int((maxvalue/minvalue)/2.)+1
@@ -650,18 +795,51 @@
cubit.cmd('set journal on')
return factor,minvalue,inv_length
+def ording_z(nodes):
+ def get_z(node):
+ x,y,z = cubit.get_nodal_coordinates(node)
+ return z
+ d = [(get_z(node), node) for node in nodes]
+ d.sort()
+ return [x[1] for x in d]
-def merge_node_4(n1,n2,n3,n4):
- factor,minvalue,inv_length=prepare_equivalence_4(n1,n2,n3,n4)
-
- for k in inv_length.keys()[:-1]:
- if len(inv_length[k]) > 0:
- cmd='equivalence node '+' '.join(' '.join(str(n) for n in x) for x in inv_length[k])+' tolerance '+str(k*factor+minvalue/2.)
+def merge_node_4(n1,n2,n3,n4,newmethod=True):
+ if newmethod:
+ print "merge node 4 side"
+ n1o=ording_z(n1)
+ n2o=ording_z(n2)
+ n3o=ording_z(n3)
+ n4o=ording_z(n4)
+ for ln in zip(n1o,n2o,n3o,n4o):
+ cmd='equivalence node '+' '.join(str(n) for n in ln) +' tolerance 10000 '
cubit.cmd(cmd)
- print 'equivalence '+str(len(inv_length[k]))+' couples of nodes - tolerance '+str(k*factor+minvalue/2.)
+ #
+ #allnodes=n1+n2+n3+n4
+ #print allnodes
+ #for n in allnodes:
+ # print n
+ #cmd='equivalence node '+' '.join(str(n) for n in allnodes) +' tolerance 10 '
+ #cubit.cmd(cmd)
+ else:
+ factor,minvalue,inv_length=prepare_equivalence_4(n1,n2,n3,n4)
+
+ for k in inv_length.keys()[:-1]:
+ if len(inv_length[k]) > 1:
+ try:
+ for x in inv_length[k]:
+ if type(x) is not list:
+ x=[x]
+ else:
+ pass
+ cmd='equivalence node '+' '.join(' '.join(str(n) for n in x) )+' tolerance '+str(k*factor+minvalue/2.)
+ except:
+ print k,"***************************************** s"
+ print inv_length[k]
+
+ cubit.cmd(cmd)
+ print 'equivalence '+str(len(inv_length[k]))+' couples of nodes - tolerance '+str(k*factor+minvalue/2.)
+ if len(inv_length[k]) == 1:
+ cmd='equivalence node '+' '.join(' '.join(str(n) for n in inv_length[k]))+' tolerance '+str(k*factor+minvalue/2.)
+ cubit.cmd(cmd)
+ print 'equivalence '+str(len(inv_length[k]))+' couples of nodes - tolerance '+str(k*factor+minvalue/2.)
-
-
-
-
-
\ No newline at end of file
Modified: seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/local_volume.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/local_volume.py 2013-09-26 00:49:25 UTC (rev 22849)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/local_volume.py 2013-09-26 00:57:45 UTC (rev 22850)
@@ -31,6 +31,145 @@
print 'error importing cubit, check if cubit is installed'
pass
+numpy = start.start_numpy()
+
+def check_orientation(grdfileNAME):
+
+ try:
+ grdfile = open(grdfileNAME, 'r')
+ print 'reading ',grdfileNAME
+ except:
+ txt='check_orintation ->error reading: '+ str( grdfile )
+ raise Exception(txt)
+ diff=1
+ txt=grdfile.readline()
+ x0,y0,z=map(float,txt.split())
+ while diff>0:
+ try:
+ txt=grdfile.readline()
+ except:
+ break
+ x,y,z=map(float,txt.split())
+ diff=x-x0
+ x0=x
+ diff=y-y0
+ if diff>0:
+ orientation= 'SOUTH2NORTH'
+ else:
+ orientation= 'NORTH2SOUTH'
+ grdfile.close()
+ return orientation
+
+def read_irregular_surf(filename):
+ "read irregular grid"
+ try:
+ xyz = numpy.loadtxt(filename)
+ except:
+ txt='error reading '+filename
+ raise Exception(txt)
+ gridpoints = xyz[:,0:2]
+ z = xyz[:,2]
+ return gridpoints,z
+
+def get_interpolated_elevation(point,gridpoints,z,k=1):
+ """for x0 and y0 return the interpolated z point of a irregular x,y,z grid
+ K=1 nearest
+ K>1 number of points used in a inverse distance weighted interpolation
+
+ point=(x0,y0)
+ gridpoints=numpy.array([[x1,y1],[x2,y2],...)
+
+ """
+ dist=numpy.sum((point-gridpoints)**2,axis=1)
+ zindex=dist.argsort()[:k]
+ kmindist=dist[zindex]
+ w=1/kmindist
+ w/=w.sum()
+ zi = numpy.dot(w.T, z[zindex])
+ return zi
+
+def create_grid(xmin,xmax,ymin,ymax,xstep,ystep):
+ """create regular grid with xmin,xmax by xstep and ymin,ymax by ystep"""
+ x,y=numpy.mgrid[xmin:xmax+xstep/2.:xstep,ymin:ymax+ystep/2.:ystep] #this includes the bounds
+ gridpoints = numpy.vstack([x.ravel(), y.ravel()]).T
+ return x,y,gridpoints
+
+
+def process_surfacefiles(iproc,nx,ny,nstep,grdfile,unit,lat_orientation):
+ from utilities import geo2utm
+ numpy = start.start_numpy()
+ elev=numpy.zeros([nx,ny],float)
+ coordx=numpy.zeros([nx,ny],float)
+ coordy=numpy.zeros([nx,ny],float)
+ icoord=0
+
+
+ lat_orientation=check_orientation(grdfile)
+
+ try:
+ grdfile = open(grdfile, 'r')
+ #print 'reading ',grdfile
+ except:
+ txt='error reading: '+ str( grdfile )
+ raise Exception(txt)
+
+
+ if lat_orientation is 'SOUTH2NORTH':
+ rangey=range(0,ny)
+ else:
+ rangey=range(ny-1,-1,-1)
+ lat_orientation='NORTH2SOUTH'
+ print lat_orientation
+ for iy in rangey:
+ for ix in range(0,nx):
+ txt=grdfile.readline()
+ try:
+ if len(txt) != 0:
+ x,y,z=map(float,txt.split())
+ if iy%nstep == 0 and ix%nstep == 0:
+ icoord=icoord+1
+ x_current,y_current=geo2utm(x,y,unit)
+ jx=min(nx-1,ix/nstep)
+ jy=min(ny-1,iy/nstep)
+ coordx[jx,jy]=x_current
+ coordy[jx,jy]=y_current
+ elev[jx,jy]=z
+ except:
+ print 'error reading point ',iy*nx+ix,txt, grdfile.name, ' proc '
+ raise NameError, 'error reading point'
+
+
+ if (nx)*(ny) != icoord:
+ if iproc == 0: print 'error in the surface file '+grdfile.name
+ if iproc == 0: print 'x points ' +str(nx)+ ' y points ' +str(ny)+ ' tot points '+str((nx)*(ny))
+ if iproc == 0: print 'points read in '+grdfile.name+': '+str(icoord)
+ raise NameError
+
+ grdfile.close()
+
+ return coordx,coordy,elev
+
+
+
+
+
+
+
+def process_irregular_surfacefiles(iproc,nx,ny,xmin,xmax,ymin,ymax,xstep,ystep,grdfile,unit_surf,lat_orientation):
+ gridpoints,z=read_irregular_surf(grdfile)
+ coordx,coordy,points=create_grid(xmin,xmax,ymin,ymax,xstep,ystep)
+
+ elev = numpy.empty([len(points)])
+ for i in xrange(len(points)):
+ elev[i] = get_interpolated_elevation(points[i],gridpoints,z,k=4)
+
+ coordx.shape=(nx,ny)
+ coordy.shape=(nx,ny)
+ elev.shape=(nx,ny)
+
+ return coordx,coordy,elev
+
+
def read_grid(filename=None):
import sys
import start as start
@@ -40,7 +179,7 @@
cfg = start.start_cfg(filename=filename)
from utilities import geo2utm
- #
+ #if cfg.irregulargridded_surf==True then cfg.nx and cfg.ny are the desired number of point along the axis....
if cfg.nx and cfg.ny:
nx=cfg.nx
ny=cfg.ny
@@ -61,182 +200,53 @@
ny= int((cfg.latitude_max-cfg.latitude_min)/ystep)+1
nstep=1
#
+
+ if cfg.irregulargridded_surf:
+ xt,xstep=numpy.linspace(cfg.xmin, cfg.xmax, num=nx, retstep=True)
+ yt,ystep=numpy.linspace(cfg.ymin, cfg.ymax, num=ny, retstep=True)
+
elev=numpy.zeros([nx,ny,cfg.nz],float)
- coordx=numpy.zeros([nx,ny],float)
- coordy=numpy.zeros([nx,ny],float)
#
if cfg.bottomflat:
elev[:,:,0] = cfg.depth_bottom
bottomsurface=1
else:
bottomsurface=0
- #
- for inz in range(bottomsurface,cfg.nz):
- try:
- grdfile = open(cfg.filename[inz-bottomsurface], 'r')
- print 'reading ',cfg.filename[inz-bottomsurface]
- except:
- txt='error reading: '+ str( cfg.filename[inz-bottomsurface] )
- raise NameError, txt
+
+ for inz in range(bottomsurface,cfg.nz-1):
+ grdfilename=cfg.filename[inz-bottomsurface]
+
+ if cfg.irregulargridded_surf:
+ coordx,coordy,elev_1=process_irregular_surfacefiles(iproc,nx,ny,cfg.xmin,cfg.xmax,cfg.ymin,cfg.ymax,xstep,ystep,grdfile)
+ else:
+ coordx,coordy,elev_1=process_surfacefiles(iproc,nx,ny,nstep,grdfilename,cfg.unit,cfg.lat_orientation)
+ elev[:,:,inz]=elev_1[:,:]
#
- icoord=0
- for iy in range(0,ny):
- for ix in range(0,nx):
- txt=grdfile.readline()
- try:
- if len(txt) != 0:
- x,y,z=map(float,txt.split())
- if iy%nstep == 0 and ix%nstep == 0:
- icoord=icoord+1
- x_current,y_current=geo2utm(x,y,cfg.unit)
- jx=min(nx-1,ix/nstep)
- jy=min(ny-1,iy/nstep)
- coordx[jx,jy]=x_current
- coordy[jx,jy]=y_current
- elev[jx,jy,inz]=z
- except:
- print 'error reading point ',iy*cfg.nx+ix,txt, cfg.filename[inz-bottomsurface], ' proc ',iproc
- raise NameError, 'error reading point'
- #
- if (nx)*(ny) != icoord:
- if iproc == 0: print 'error in the surface file '+cfg.filename[inz-bottomsurface]
- if iproc == 0: print 'x points ' +str(nx)+ ' y points ' +str(ny)+ ' tot points '+str((nx)*(ny))
- if iproc == 0: print 'points read in '+cfg.filename[inz-bottomsurface]+': '+str(icoord)
- raise NameError
-
- #if iproc == 0: print 'end of reading grd ascii file '+cfg.filename[inz-bottomsurface]+' '+str(icoord)+ ' points'
- grdfile.close()
-
+ inz=cfg.nz-1 #last surface
+ if cfg.sea:
+ elev[:,:,inz]=elev[:,:,inz-1]
+ else:
+ #try:
+ grdfile = cfg.filename[inz-bottomsurface]
+ print 'reading ',cfg.filename[inz-bottomsurface]
+ if cfg.irregulargridded_surf:
+ coordx,coordy,elev_1=process_irregular_surfacefiles(iproc,nx,ny,cfg.xmin,cfg.xmax,cfg.ymin,cfg.ymax,xstep,ystep,grdfile)
+ else:
+ coordx,coordy,elev_1=process_surfacefiles(iproc,nx,ny,nstep,grdfile,cfg.unit,cfg.lat_orientation)
+ elev[:,:,inz]=elev_1[:,:]
+ #except:
+ # txt='error reading: '+ str( cfg.filename[inz-bottomsurface] )
+ # raise NameError, txt
+
+
+ if cfg.subduction:
+ print 'subduction'
+ top=elev[:,:,inz]
+ slab=elev[:,:,inz-1]
+ subcrit=numpy.abs(top-slab)<cfg.subduction_thres
+ top[subcrit]=slab[subcrit]+cfg.subduction_thres
+ print len(top[subcrit])
+ elev[:,:,inz]=top
return coordx,coordy,elev,nx,ny
-
-def extract_volume(xmin,ymin,xmax,ymax,coordx,coordy,elev,nx,ny,filename=None):
- import sys
- import start as start
- #
- mpiflag,iproc,numproc,mpi = start.start_mpi()
- #
- numpy = start.start_numpy()
- cfg = start.start_cfg(filename=filename)
-
- from utilities import geo2utm
- #
- rxstep=coordx[1,0]-coordx[0,0]
- rystep=coordy[0,1]-coordy[0,0]
-
- nxmin_cpu=min(0,int((x0-cfg.xmin)/rxstep)+1-10)
- nymin_cpu=min(0,int((y0-cfg.ymin)/rxstep)+1-10)
- nxmax_cpu=min(nx,int((x0-cfg.xmin)/rystep)+1+10)
- nymax_cpu=min(ny,int((y0-cfg.ymin)/rystep)+1+10)
- #
- #
- icurve=0
- isurf=0
- ivertex=0
- #
- #create vertex
- last_surface=cubit.get_last_id('surface')
- for inz in range(0,cfg.nz):
- if cfg.bottomflat and inz == 0: #bottom layer
-
- x_current,y_current=geo2utm(coordx[nxmin_cpu,nymin_cpu],coordy[nxmin_cpu,nymin_cpu],cfg.unit)
- cubitcommand= 'create vertex '+ str( x_current )+ ' ' + str( y_current) +' '+ str( cfg.depth_bottom )
- cubit.cmd(cubitcommand)
- #
- x_current,y_current=geo2utm(coordx[nxmin_cpu,nymax_cpu],coordy[nxmin_cpu,nymax_cpu],cfg.unit)
- cubitcommand= 'create vertex '+ str( x_current )+ ' ' + str( y_current) +' '+ str( cfg.depth_bottom )
- cubit.cmd(cubitcommand)
- #
- x_current,y_current=geo2utm(coordx[nxmax_cpu,nymax_cpu],coordy[nxmax_cpu,nymax_cpu],cfg.unit)
- cubitcommand= 'create vertex '+ str( x_current )+ ' ' + str( y_current) +' '+ str( cfg.depth_bottom )
- cubit.cmd(cubitcommand)
- #
- x_current,y_current=geo2utm(coordx[nxmax_cpu,nymin_cpu],coordy[nxmax_cpu,nymin_cpu],cfg.unit)
- cubitcommand= 'create vertex '+ str( x_current )+ ' ' + str( y_current) +' '+ str( cfg.depth_bottom )
- cubit.cmd(cubitcommand)
- #
- cubitcommand= 'create surface vertex 1 2 3 4'
- cubit.cmd(cubitcommand)
- #
- isurf = isurf + 1
-
- else:
- vertex=[]
-
- for iy in range(nymin_cpu,nymax_cpu+1):
- ivx=0
- for ix in range(nxmin_cpu,nxmax_cpu+1):
- zvertex=elev[ix,iy,inz]
- x_current,y_current=geo2utm(coordx[ix,iy],coordy[ix,iy],cfg.unit)
- #
- vertex.append(' Position '+ str( x_current ) +' '+ str( y_current )+' '+ str( zvertex ) )
- #
- print iproc, 'vertex created....'
- n=max(nx,ny)
- uline=[]
- vline=[]
- iv=0
-
- cubit.cmd("set info off")
- cubit.cmd("set echo off")
- cubit.cmd("set journal off")
-
- for iy in range(0,nymax_cpu-nymin_cpu+1):
- positionx=''
- for ix in range(0,nxmax_cpu-nxmin_cpu+1):
- positionx=positionx+vertex[iv]
- iv=iv+1
- command='create curve spline '+positionx
- cubit.cmd(command)
- uline.append( cubit.get_last_id("curve") )
- for ix in range(0,nxmax_cpu-nxmin_cpu+1):
- positiony=''
- for iy in range(0,nymax_cpu-nymin_cpu+1):
- positiony=positiony+vertex[ix+iy*(nxmax_cpu-nxmin_cpu+1)]
- command='create curve spline '+positiony
- cubit.cmd(command)
- vline.append( cubit.get_last_id("curve") )
- #
- cubit.cmd("set info "+cfg.cubit_info)
- cubit.cmd("set echo "+cfg.echo_info)
- cubit.cmd("set journal "+cfg.jou_info)
- #
- #
- print iproc,'line created....'
- umax=max(uline)
- umin=min(uline)
- vmax=max(vline)
- vmin=min(vline)
- cubitcommand= 'create surface net u curve '+ str( umin )+' to '+str( umax )+ ' v curve '+ str( vmin )+ ' to '+str( vmax )+' heal'
- cubit.cmd(cubitcommand)
- command = "del curve all"
- cubit.cmd(command)
- isurf=isurf+1
- #
- #
- cubitcommand= 'del vertex all'
- cubit.cmd(cubitcommand)
- #cubit_error_stop(iproc,cubitcommand,ner)
- cubitcommand= 'del curve all'
- cubit.cmd(cubitcommand)
- #
- last_surface_2=cubit.get_last_id('surface')
- #
- for inz in range(1,cfg.nz):
- #!cubit cmd
- cubitcommand= 'create volume loft surface '+ str( inz+1 )+' '+str( inz )
- cubit.cmd(cubitcommand)
- #cubit_error_stop(iproc,cubitcommand,ner)
- isurf=isurf+6
- cubitcommand= 'del surface '+str(last_surface+1)+' to '+ str( last_surface_2 )
- cubit.cmd(cubitcommand)
- #cubit_error_stop(iproc,cubitcommand,ner)
- #
- #
- cubit.cmd("set info "+cfg.cubit_info)
- cubit.cmd("set echo "+cfg.echo_info)
- cubit.cmd("set journal "+cfg.jou_info)
- command = "compress all"
- cubit.cmd(command)
-
\ No newline at end of file
Modified: seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/menu.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/menu.py 2013-09-26 00:49:25 UTC (rev 22849)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/menu.py 2013-09-26 00:57:45 UTC (rev 22850)
@@ -68,37 +68,30 @@
collect some cubit files and merge in a single free mesh cubitfile
GEOCUBIT.py --collect --merge --meshfiles=[list of files] --cpux=N --cpuy=N (--rangecpux=[cpuxmin,cpuxmax], --rangecpuy=[cpuymin,cpuymax])
+
+ collect some cubit files and merge in a single free mesh cubitfile (and decimate!!! (refine by 2))
+ GEOCUBIT.py --collect --decimate --merge --meshfiles=[list of files] --cpux=N --cpuy=N (--rangecpux=[cpuxmin,cpuxmax], --rangecpuy=[cpuymin,cpuymax])
+
+
collect a single free mesh cubitfile and refine the hex inside some curve (ex. basin)
GEOCUBIT.py --collect --meshfiles=[list of files] --curverefining=[list of SAT files]
- export a cubit mesh file (with blocks defined following the note) in a SPECFEM3D_Cartesian mesh
+ export a cubit mesh file (with blocks defined following the note) in a SPECFEM3D_SESAME mesh
GEOCUBIT.py --export2SPECFEM3D --meshfiles=[filename] (--listblock=block1,block2,..,blockN --listflag=[list of specfem flag, i.e. --listflag=1,2,3,-1])
"""
try:
- opts, args = getopt.getopt(sys.argv[1:], "sjmohbp1", ["SEMoutput=","qlog","mfast","curverefining=","output=","rangecpux=","rangecpuy=","equivalence","listflag=","listblock=","cpux=","cpuy=","exofiles=","partitioner","plane","x1=","x2=","x3=","x4=","unit=","chkcfg","mat=","merge_tolerance=","export2SPECFEM3D=","mesh","chklib","cfg=","job=","basin","help", "id_proc=", "surface=","script","jou","strat","MPI","regulargrid=",'skin=',"build_surface","build_volume","merge1","merge2","merge","collect","meshfiles="])
+ opts, args = getopt.getopt(sys.argv[1:], "sjmohbp1", ["hex27","cpml_size=","top_absorbing","cpml","decimate","addsea","SEMoutput=","qlog","mfast","curverefining=","output=","rangecpux=","rangecpuy=","equivalence","listflag=","listblock=","cpux=","cpuy=","exofiles=","partitioner","plane","x1=","x2=","x3=","x4=","unit=","chkcfg","mat=","merge_tolerance=","export2SPECFEM3D","mesh","chklib","cfg=","job=","basin","help", "id_proc=", "surface=","script","jou","strat","MPI","regulargrid=",'skin=',"build_surface","build_volume","merge1","merge2","merge","collect","meshfiles="])
print opts, args
-except getopt.GetoptError,errmsg:
- if str(errmsg) == 'option --export2SPECFEM3D requires argument':
- for i,xop in enumerate(sys.argv[1:]):
- if 'export2SPECFEM3D' in xop:
- sys.argv[i+1]=sys.argv[i+1]+'=.'
- try:
- opts, args = getopt.getopt(sys.argv[1:], "sjmohbp1", ["SEMoutput=","qlog","mfast","curverefining=","output=","rangecpux=","rangecpuy=","equivalence","listflag=","listblock=","cpux=","cpuy=","exofiles=","partitioner","plane","x1=","x2=","x3=","x4=","unit=","chkcfg","mat=","merge_tolerance=","export2SPECFEM3D=","mesh","chklib","cfg=","job=","basin","help", "id_proc=", "surface=","script","jou","strat","MPI","regulargrid=",'skin=',"build_surface","build_volume","merge1","merge2","merge","collect","meshfiles="])
- print opts, args
- except getopt.GetoptError,errmsg:
- print str(errmsg)
- usage()
- sys.exit(2)
- else:
- print str(errmsg)
- usage()
- sys.exit()
-except AttributeError:
- opts=None
- args=None
- print opts, args
+except Exception,e:
+ #if 'argv' in e:
+ # pass
+ #else:
+ # if '__console__' != __name__: #if you are in the cubitGUI/pythonprompt you don't want usage()
+ # usage()
+ # print e
+ sys.exit()
output='totalmesh_merged'
SPECFEM3D_output_dir='.'
@@ -136,17 +129,35 @@
cpuxmax=None
cpuymax=None
curverefining=False
+add_sea=False
+decimate=False
+cpml=False
+cpml_size=False
+top_absorbing=False
-
-
-
qlog=False
+hex27=False
if opts:
for o, value in opts:
#print o,value
+ if o in ('--hex27'):
+ hex27=True
+ if o in ('--cpml'):
+ cpml=True
+ if '--cpml_size' in o:
+ cpml=True
+ cpml_size=float(value)
+ if o in ('--top_absorbing'):
+ cpml=True
+ top_absorbing=True
+ if '--cpml_size' in o:
+ cpml=True
+ cpml_size=float(value)
+ if o in ('--decimate'):
+ decimate=True
if o in ('--partitioner'):
create_partitioner=True
if o == ('--surface'):
@@ -209,6 +220,13 @@
sys.exit()
if o in ("--cfg"):
cfg_name=value
+ try:
+ if open(cfg_name): pass
+ except IOError,e:
+ print 'error opening ',cfg_name
+ print e
+ import sys
+ sys.exit()
if o == ('--surface'):
surface=True
surface_name=value
@@ -231,12 +249,6 @@
exofiles=value
if o in ("--export2SPECFEM3D"):
export2SPECFEM3D=True
- SPECFEM3D_output_dir=value
- import os
- try:
- os.makedirs(SPECFEM3D_output_dir)
- except OSError:
- pass
if o in ("--merge_tolerance") and o != '--merge' and o != '--merge2' and o != '--merge1':
merge_tolerance=map(float,value.split(','))
if o in ("--mat"):
@@ -265,8 +277,15 @@
output=value
if o in ("--curverefining"):
curverefining=value.split(',')
- if o in ("SEMoutput"):
+ if o in ("--SEMoutput"):
SPECFEM3D_output_dir=value
+ import os
+ try:
+ os.makedirs(SPECFEM3D_output_dir)
+ except OSError:
+ pass
+ if o in ("--addsea"):
+ add_sea=True
print cpuxmax,cpuymax
if cpuymax:
pass
@@ -281,6 +300,18 @@
else:
cpuxmax=1
print cpuxmax,cpuymax
+
+ if cpml:
+ if not cpml_size:
+ print 'specify the size of the cpml boundaries'
+ import sys
+ sys.exit()
+ elif cpml_size<=0:
+ print 'cpml size negative, please check the parameters'
+ import sys
+ sys.exit()
+
+
if chkcfg==True:
import start as start
cfg=start.start_cfg()
@@ -312,4 +343,5 @@
elif opts == []:
print __name__
usage()
- raise AttributeError('no options')
+ import sys
+ sys.exit()
Modified: seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/mesh_volume.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/mesh_volume.py 2013-09-26 00:49:25 UTC (rev 22849)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/mesh_volume.py 2013-09-26 00:57:45 UTC (rev 22850)
@@ -41,16 +41,23 @@
#
if cfg.map_meshing_type == 'regularmap':
mesh_layercake_regularmap(filename=filename)
- #elif cfg.map_meshing_type == 'partitioner':
- # mesh_partitioner()
else:
print 'error: map_meshing_type ', cfg.map_meshing_type,' not implemented'
+### AAA edited 8/28
def mesh_layercake_regularmap(filename=None):
import sys,os
import start as start
+ mpiflag,iproc,numproc,mpi = start.start_mpi()
+ from utilities import importgeometry,savemesh,get_v_h_list,cubit_command_check
#
+ numpy = start.start_numpy()
+ cfg = start.start_cfg(filename=filename)
+ from math import sqrt
+ from sets import Set
+
+ #
class cubitvolume:
def __init__(self,ID,intervalv,centerpoint,dimension):
self.ID=ID
@@ -66,14 +73,7 @@
return cmp(x.centerpoint,y.centerpoint)
#
#
- mpiflag,iproc,numproc,mpi = start.start_mpi()
- from utilities import importgeometry,savemesh,get_v_h_list
#
- numpy = start.start_numpy()
- cfg = start.start_cfg(filename=filename)
- from math import sqrt
- from sets import Set
- #
list_vol=cubit.parse_cubit_list("volume","all")
if len(list_vol) != 0:
pass
@@ -83,7 +83,7 @@
#
command = 'composite create curve all'
cubit.cmd(command)
- print 'NO CRITICAL ERROR: "No valid composites can be created from the specified curves." is not critical. \n It means that your model is clean and you don"t need a virtual geometry'
+ print '###"No valid composites can be created from the specified curves." is NOT a critical ERROR.'
#
command = "compress all"
cubit.cmd(command)
@@ -113,7 +113,8 @@
#
#
#interval assignement
- surf_or,surf_vertical,list_curve_or,list_curve_vertical,bottom,top = get_v_h_list(list_vol)
+ surf_or,surf_vertical,list_curve_or,list_curve_vertical,bottom,top = get_v_h_list(list_vol,chktop=cfg.chktop)
+ print 'vertical surfaces: ',surf_vertical
for k in surf_vertical:
command = "surface "+str(k)+" scheme submap"
@@ -121,16 +122,47 @@
for k in surf_or:
command = "surface "+str(k)+" scheme "+cfg.or_mesh_scheme
cubit.cmd(command)
- for k in list_curve_or:
+ #
+ ucurve,vcurve=get_uv_curve(list_curve_or)
+ schemepave=False
+ #
+ ucurve_interval={}
+ for k in ucurve:
length=cubit.get_curve_length(k)
interval=int(2*round(.5*length/cfg.size,0))
+ ucurve_interval[k]=interval
command = "curve "+str(k)+" interval "+str(interval)
cubit.cmd(command)
#cubit_error_stop(iproc,command,ner)
command = "curve "+str(k)+" scheme equal"
cubit.cmd(command)
#cubit_error_stop(iproc,command,ner)
- #
+ if max(ucurve_interval.values()) != min(ucurve_interval.values()):
+ schemepave=True
+ print 'mesh scheme is set to pave'
+ for sk in surf_or:
+ command = "surface "+str(sk)+" scheme pave"
+ cubit.cmd(command)
+ #
+ vcurve_interval={}
+ for k in vcurve:
+ length=cubit.get_curve_length(k)
+ interval=int(2*round(.5*length/cfg.size,0))
+ vcurve_interval[k]=interval
+ command = "curve "+str(k)+" interval "+str(interval)
+ cubit.cmd(command)
+ #cubit_error_stop(iproc,command,ner)
+ command = "curve "+str(k)+" scheme equal"
+ cubit.cmd(command)
+ #cubit_error_stop(iproc,command,ner)
+
+ if max(vcurve_interval.values()) != min(vcurve_interval.values()):
+ print 'mesh scheme is set to pave'
+ schemepave=True
+ for sk in surf_or:
+ command = "surface "+str(sk)+" scheme pave"
+ cubit.cmd(command)
+ #
for s in surf_vertical:
lcurve=cubit.get_relatives("surface",s,"curve")
interval_store=[]
@@ -157,18 +189,30 @@
command = "curve "+' '.join(str(iv[0]) for iv in interval_store)+" scheme equal"
cubit.cmd(command)
#cubit_error_stop(iproc,command,ner)
+ command = "surface "+str(s)+" scheme submap"
+ cubit.cmd(command)
+
#cubit_error_stop(iproc,command,ner)
#
#meshing
- if cfg.or_mesh_scheme == 'pave':
+ if cfg.or_mesh_scheme == 'pave' or schemepave:
command='mesh surf '+' '.join(str(t) for t in top)
- cubit.cmd(command)
+ status=cubit_command_check(iproc,command,stop=True)
+ #cubit.cmd(command)
elif cfg.or_mesh_scheme == 'map':
command='mesh surf '+' '.join(str(t) for t in bottom)
- cubit.cmd(command)
+ status=cubit_command_check(iproc,command,stop=True)
+ #cubit.cmd(command)
for id_volume in range(nvol-1,-1,-1):
command = "mesh vol "+str(vol[id_volume].ID)
- cubit.cmd(command)
+ status=cubit_command_check(iproc,command,stop=False)
+ if not status:
+ for s in surf_vertical:
+ command_surf="mesh surf "+str(s)
+ cubit.cmd(command_surf)
+ command_set_meshvol='volume all redistribute nodes on\nvolume all autosmooth target off\nvolume all scheme Sweep Vector 0 0 -1\nvolume all sweep smooth Auto\n'
+ status=cubit_command_check(iproc,command_set_meshvol,stop=False)
+ status=cubit_command_check(iproc,command,stop=True)
#
#smoothing
@@ -216,7 +260,6 @@
cubitcommand= 'del mesh vol '+str(vol[-1].ID)+ ' propagate'
cubit.cmd(cubitcommand)
s1=Set(list_curve_vertical)
- print s1
command = "group 'list_curve_tmp' add curve "+"in vol "+str(vol[-1].ID)
cubit.cmd(command)
group=cubit.get_id_from_name("list_curve_tmp")
@@ -224,9 +267,7 @@
command = "delete group "+ str(group)
cubit.cmd(command)
s2=Set(list_curve_tmp)
- print s2
lc=list(s1 & s2)
- print lc
#
cubitcommand= 'curve '+' '.join(str(x) for x in lc)+' interval '+str(cfg.actual_vertical_interval_top_layer)
cubit.cmd(cubitcommand)
@@ -245,7 +286,7 @@
import boundary_definition
entities=['face']
print iproc, 'hex block definition...'
- boundary_definition.define_bc(entities,parallel=True,cpux=cfg.cpux,cpuy=cfg.cpuy,cpuxmin=0,cpuymin=0)
+ boundary_definition.define_bc(entities,parallel=True,cpux=cfg.cpux,cpuy=cfg.cpuy,cpuxmin=0,cpuymin=0,optionsea=False)
#save mesh
print iproc, 'untangling...'
@@ -261,6 +302,7 @@
def refinement(nvol,vol,filename=None):
import start as start
cfg = start.start_cfg(filename=filename)
+ from utilities import get_v_h_list
#
#vertical refinement
#for nvol = 3
@@ -280,14 +322,17 @@
if cfg.ntripl != 0:
if len(cfg.refinement_depth) != 0:
#get the topo surface....
- surf=cubit.get_relatives('volume',vol[nvol+1-2].ID,'surface')
- zstore=[-1,-999999999]
- for s in surf:
- c=cubit.get_center_point('surface',s)
- z=c[2]
- if z > zstore[1]:
- zstore=[s,z]
- tsurf=zstore[0]
+ #surf=cubit.get_relatives('volume',vol[nvol-1].ID,'surface')
+ #zstore=[-1,-999999999]
+ #for s in surf:
+ # c=cubit.get_center_point('surface',s)
+ # z=c[2]
+ # print s,z
+ # if z > zstore[1]:
+ # zstore=[s,z]
+ #tsurf=zstore[0]
+ _,_,_,_,_,tsurf = get_v_h_list([vol[nvol-1].ID])
+ tsurf=' '.join(str(x) for x in tsurf)
for idepth in cfg.refinement_depth:
cubitcommand= 'refine node in surf '+str(tsurf)+' numsplit 1 bias 1.0 depth '+str(idepth)
cubit.cmd(cubitcommand)
@@ -313,18 +358,18 @@
cubitcommand= 'refine hex in vol '+txt
else:
#refinement on the top surface
-
- surf=cubit.get_relatives('volume',vol[ir-2].ID,'surface')
- zstore=[-1,-999999999]
- for s in surf:
- c=cubit.get_center_point('surface',s)
- z=c[2]
- if z > zstore[1]:
- zstore=[s,z]
+ _,_,_,_,_,tsurf = get_v_h_list([vol[ir-2].ID])
+ tsurf=' '.join(str(x) for x in tsurf)
idepth=1
- cubitcommand= 'refine node in surf '+str(zstore[0])+' numsplit 1 bias 1.0 depth '+str(idepth)
+ cubitcommand= 'refine node in surf '+str(tsurf)+' numsplit 1 bias 1.0 depth '+str(idepth)
cubit.cmd(cubitcommand)
+ if not nvol and cfg.volume_type == 'verticalsandwich_volume_ascii_regulargrid_mpiregularmap':
+ # AAA
+ # Volume 2 is always in between the 2nd and 3rd vertical surfaces from the left
+ cubitcommand = "refine node in volume 2 numsplit 1 depth 0"
+ cubit.cmd(cubitcommand)
+ # END AAA
@@ -441,3 +486,36 @@
import sys
#sys.exit()
+def get_uv_curve(list_curve_or):
+ import math
+ import numpy as np
+ klen={}
+ for curve in list_curve_or:
+ vertex_list = cubit.get_relatives("curve", curve, "vertex")
+ coord0=cubit.get_center_point('vertex', vertex_list[0])
+ coord1=cubit.get_center_point('vertex', vertex_list[1])
+ klen[curve]=np.array(coord1)-np.array(coord0)
+ #
+ l0=list_curve_or[0]
+ c0=klen[l0]
+ angles={}
+ angles[l0]=0
+ for curve in list_curve_or[1:]:
+ c1=klen[curve]
+ angletmp=np.dot(c0,c1)/(np.dot(c0,c0)**.5*np.dot(c1,c1)**.5)
+ if -1 < angletmp < 1:
+ angle=math.sin(np.arccos(angletmp))
+ else:
+ angle=0.
+ angles[curve]=angle
+ a=angles.values()
+ diff=max(a)-min(a)
+ ucurve=[]
+ vcurve=[]
+ for curve in list_curve_or:
+ if -diff < angles[curve] < diff:
+ ucurve.append(curve)
+ else:
+ vcurve.append(curve)
+ #
+ return ucurve,vcurve
\ No newline at end of file
Modified: seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/read_parameter_cfg.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/read_parameter_cfg.py 2013-09-26 00:49:25 UTC (rev 22849)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/read_parameter_cfg.py 2013-09-26 00:57:45 UTC (rev 22850)
@@ -46,7 +46,7 @@
import sys
sys.exit()
#
- from utilities import geo2utm #here I can use pyproj but I prefere to include a function in pure python in order to avoid an additional installation
+ from utilities import geo2utm,get_cubit_version #here I can use pyproj but I prefere to include a function in pure python in order to avoid an additional installation
#
#
import ConfigParser
@@ -119,6 +119,7 @@
#CONSTANTS
dcfg['osystem'] = 'linux'
dcfg['debug_cfg']=False
+ dcfg['version_cubit']=get_cubit_version()
dcfg['checkbound']=False
dcfg['top_partitioner'] = 10000
dcfg['tres']=0.3 #if n is the vertical component of the normal at a surface pointing horizontally, when -tres < n < tres then the surface is vertical
@@ -157,10 +158,21 @@
dcfg['filename']=None
dcfg['actual_vertical_interval_top_layer']=1
dcfg['coarsening_top_layer']=False
+ dcfg['refineinsidevol']=False
+ dcfg['sea']=False
+ dcfg['seaup']=False
+ dcfg['sea_level']=False
+ dcfg['sea_threshold']=False
+ dcfg['subduction']=True
+ dcfg['subduction_thres']=500
+ dcfg['debugsurface']=False #if true it creates only the surface not the lofted volumes
+ dcfg['lat_orientation']=False
+ dcfg['irregulargridded_surf']=False
+ dcfg['chktop']=False
+ dcfg['smoothing']=False
+ dcfg['ntripl']=0
-
-
dcfg['nsurf'] = None
if cfgname:
config.read(cfgname)
@@ -173,7 +185,7 @@
dcfg.update(d)
except:
pass
-
+
if dcfg['nsurf']:
surface_name=[]
num_x=[]
@@ -242,6 +254,14 @@
except:
pass
+ if dcfg['sea']:
+ if not dcfg['sea_level']: dcfg['sea_level']=0
+ if not dcfg['sea_threshold']: dcfg['sea_threshold']=-200
+ dcfg['actual_vertical_interval_top_layer']=1
+ dcfg['coarsening_top_layer']=True
+
+
+ dcfg['optionsea']={'sea':dcfg['sea'],'seaup':dcfg['seaup'],'sealevel':dcfg['sea_level'],'seathres':dcfg['sea_threshold']}
cfg=attrdict(dcfg)
if menu:
@@ -312,7 +332,16 @@
except:
pass
+ try:
+ if isinstance(cfg.iv_interval,int): cfg.iv_interval=[cfg.iv_interval]
+ except:
+ pass
+ try:
+ if isinstance(cfg.refinement_depth,int): cfg.refinement_depth=[cfg.refinement_depth]
+ except:
+ pass
+
return cfg
Deleted: seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/save_fault_nodes_elements.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/save_fault_nodes_elements.py 2013-09-26 00:49:25 UTC (rev 22849)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/save_fault_nodes_elements.py 2013-09-26 00:57:45 UTC (rev 22850)
@@ -1,84 +0,0 @@
-#!python
-# Opening fault cracks
-# Input : surface up and down.
-import cubit
-
-class fault_input:
- def __init__(self,id,surface_u,surface_d):
- self.id = id
- self.surface_u = surface_u
- self.surface_d = surface_d
- self.name = 'MESH/fault_file_'+str(id)+'.dat'
-
- quads_Aup,quads_Adp = save_cracks(self.name,self.surface_u,self.surface_d)
- #Unpacking list.
- quads_Au=unpack_list(quads_Aup)
- quads_Ad=unpack_list(quads_Adp)
-
- print 'len(Au):',len(quads_Au)
- print 'len(Ad):',len(quads_Ad)
-
- if not (len(quads_Au)==len(quads_Ad)):
- print 'Number of elements for each fauld side up and down do not concide'
- sys.exit('goodbye')
-
- save_elements_nodes(self.name,quads_Au,quads_Ad)
-
-
-def save_cracks(name,list_surface_up,list_surface_down):
- quads_fault_up = []
- quads_fault_down = []
- for surface in list_surface_up :
- quads_fault = cubit.get_surface_quads(surface)
- quads_fault_up.append(quads_fault)
- for surface in list_surface_down :
- quads_fault = cubit.get_surface_quads(surface)
- quads_fault_down.append(quads_fault)
- # TO DO : stop python properly in case fault nodes at both sides
- # do not match.
- # if len(quads_fault_u) != len(quads_fault_d): stop
- #
- # SAVING FAULT ELEMENTS AND NODES
- return quads_fault_up,quads_fault_down
-
-def unpack_list(fault_list):
- list_fault = []
- for i in range(0,len(fault_list)):
- el=list(fault_list[i])
- for j in el:
- list_fault.append(j)
- return list_fault
-
-def save_elements_nodes(name,quads_fault_u,quads_fault_d):
- fault_file = open(name,'w')
- txt =''
- list_hex=cubit.parse_cubit_list('hex','all')
- txt='%10i %10i\n' % (len(quads_fault_u),len(quads_fault_d))
- fault_file.write(txt)
-
- dic_quads_fault_u = dict(zip(quads_fault_u,quads_fault_u))
- dic_quads_fault_d = dict(zip(quads_fault_d,quads_fault_d))
-
- # FAULT SIDE DOWN
- for h in list_hex:
- faces = cubit.get_sub_elements('hex',h,2)
- for f in faces:
- if dic_quads_fault_d.has_key(f):
- nodes=cubit.get_connectivity('Face',f)
- # print 'h,fault nodes side down :',h,nodes[0],nodes[1],nodes[2],nodes[3]
- txt='%10i %10i %10i %10i %10i\n' % (h,nodes[0],\
- nodes[1],nodes[2],nodes[3])
- fault_file.write(txt)
-
- # FAULT SIDE UP
- for h in list_hex:
- faces = cubit.get_sub_elements('hex',h,2)
- for f in faces:
- if dic_quads_fault_u.has_key(f):
- nodes=cubit.get_connectivity('Face',f)
- # print 'h,fault nodes side up :',h,nodes[0],nodes[1],nodes[2],nodes[3]
- txt='%10i %10i %10i %10i %10i\n' % (h,nodes[0],\
- nodes[1],nodes[2],nodes[3])
- fault_file.write(txt)
-
- fault_file.close()
Added: seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/sea.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/sea.py (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/sea.py 2013-09-26 00:57:45 UTC (rev 22850)
@@ -0,0 +1,40 @@
+#############################################################################
+# sea.py
+# this file is part of GEOCUBIT #
+# #
+# Created by Emanuele Casarotti #
+# Copyright (c) 2008 Istituto Nazionale di Geofisica e Vulcanologia #
+# #
+#############################################################################
+# #
+# GEOCUBIT is free software: you can redistribute it and/or modify #
+# it under the terms of the GNU General Public License as published by #
+# the Free Software Foundation, either version 3 of the License, or #
+# (at your option) any later version. #
+# #
+# GEOCUBIT is distributed in the hope that it will be useful, #
+# but WITHOUT ANY WARRANTY; without even the implied warranty of #
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
+# GNU General Public License for more details. #
+# #
+# You should have received a copy of the GNU General Public License #
+# along with GEOCUBIT. If not, see <http://www.gnu.org/licenses/>. #
+# #
+#############################################################################
+
+def adjust_sea_layers(zvertex,sealevel,bathymetry,cfg):
+ if cfg.seaup:
+ if sealevel:
+ vertex=max(cfg.sea_level,vertex-cfg.sea_threshold)
+ elif bathymetry:
+ #if zvertex > cfg.sea_threshold: zvertex=max(zvertex,cfg.sea_level)+cfg.sea_threshold #move node below the topography of sea_threshold
+ vertex=vertex
+ else:
+ if sealevel and zvertex < cfg.sea_level:
+ zvertex=cfg.sea_level
+ elif bathymetry:
+ if zvertex > cfg.sea_threshold: zvertex=max(zvertex,cfg.sea_level)+cfg.sea_threshold #move node below the topography of sea_threshold
+ return zvertex
+
+
+
Modified: seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/start.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/start.py 2013-09-26 00:49:25 UTC (rev 22849)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/start.py 2013-09-26 00:57:45 UTC (rev 22850)
@@ -85,6 +85,7 @@
except:
try:
import cubit
+ import utilities
cubit.init([""])
except:
print 'error importing cubit'
@@ -122,6 +123,14 @@
cubit.cmd('comment "'+txt+'"')
cubit.cmd("set echo "+cfg.echo_info)
cubit.cmd("set info "+cfg.cubit_info)
+ version_cubit=utilities.get_cubit_version()
+ if version_cubit <= 13:
+ print 'VERSION CUBIT ',version_cubit
+ elif version_cubit > 13:
+ print 'CAVEAT:'
+ print 'VERSION CUBIT ',version_cubit
+ print 'VERSIONs of CUBIT > 13 have bugs with merge node commands and equivalence'
+ print 'the merge option is not operative with this version, please download CUBIT 13'
except:
print 'error start cubit'
sys.exit()
Modified: seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/utilities.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/utilities.py 2013-09-26 00:49:25 UTC (rev 22849)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/utilities.py 2013-09-26 00:57:45 UTC (rev 22850)
@@ -32,8 +32,15 @@
pass
+def get_cubit_version():
+ v=cubit.get_version()
+ try:
+ v=float(v[0:4])
+ except:
+ v=float(v[0:2])
+ return v
+
-
def snapshot(name=None,i=0,viewnumber=1):
"""
it takes a snapshot of the figure, following the predefined view position.
@@ -55,7 +62,38 @@
cubit.cmd(command)
return i
+def cubit_command_check(iproc,command,stop=True):
+ """
+ Run a cubit command, checking if it performs correctly.
+ If the command fails, it writes the result on a file "error_[processor number]" and stop the meshing process
+
+ iproc = process number
+ command = cubit command
+ stop = if command fails, stop the meshing process (Default: True)
+
+ return status variable (0 ok, -1 fail)
+
+ """
+ er=cubit.get_error_count()
+ cubit.cmd(command)
+ ner=cubit.get_error_count()
+ flag=0
+ if ner > er:
+ text='"Proc: '+str(iproc)+' ERROR '+str(command)+' number of error '+str(er)+'/'+str(ner)+'"'
+ cubitcommand = 'comment '+text
+ cubit.cmd(cubitcommand)
+ f=open('error_'+str(iproc)+'.log','a')
+ f.write("CUBIT ERROR: \n"+text)
+ f.close()
+ if stop: raise Exception("CUBIT ERROR: "+text)
+ flag=-1
+ return flag
+
+
+
+
def cubit_error_stop(iproc,command,ner):
+ """obsolete"""
er=cubit.get_error_count()
if er > ner:
text='"Proc: '+str(iproc)+' ERROR '+str(command)+' number of error '+str(er)+'/'+str(ner)+'"'
@@ -64,6 +102,7 @@
raise NameError, text
def cubit_error_continue(iproc,command,n_er):
+ """obsolete"""
er=cubit.get_error_count()
if er >= n_er:
text='"Proc: '+str(iproc)+' ERROR continue '+str(command)+' number of error '+str(er)+' '+str(n_er)+'"'
@@ -359,14 +398,14 @@
f.close()
-def get_v_h_list(vol_id_list):
+def get_v_h_list(vol_id_list,chktop=False):
"""return the lists of the cubit ID of vertical/horizontal surface and vertical/horizontal curves
where v/h is defined by the distance of the z normal component from the axis direction
the parameter cfg.tres is the threshold as for example if
- normal[2] >= -tres and normal[2] <= tres
+ -tres <= normal[2] <= tres
then the surface is vertical
#
- usage: surf_or,surf_vertical,list_curve_or,list_curve_vertical,bottom,top = get_v_h_list(list_vol)
+ usage: surf_or,surf_vertical,list_curve_or,list_curve_vertical,bottom,top = get_v_h_list(list_vol,chktop=False)
"""
#
tres=0.3
@@ -387,7 +426,7 @@
for k in lsurf:
normal=cubit.get_surface_normal(k)
center_point = cubit.get_center_point("surface", k)
- if normal[2] >= -1*tres and normal[2] <= tres:
+ if -1*tres <= normal[2] <= tres:
surf_vertical.append(k)
lcurve=cubit.get_relatives("surface",k,"curve")
list_curve_vertical=list_curve_vertical+list(lcurve)
@@ -400,6 +439,8 @@
list_curve_vertical.remove(x)
except:
pass
+
+ #find the top and the bottom surfaces
k=surf_or[0]
center_point = cubit.get_center_point("surface", k)[2]
center_point_top=center_point
@@ -414,14 +455,49 @@
elif center_point < center_point_bottom:
center_point_bottom=center_point
bottom=k
- surftop=list(cubit.get_adjacent_surfaces("surface", top))
+ #check that a top surface exists
+ #it assume that the z coord of the center point
+ if chktop:
+ k=lsurf[0]
+ vertical_centerpoint_top = cubit.get_center_point("surface", k)[2]
+ vertical_zmax_box_top=cubit.get_bounding_box('surface',k)[7]
+ normal_top=cubit.get_surface_normal(k)
+ top=k
+ for k in lsurf:
+ vertical_centerpoint = cubit.get_center_point("surface", k)[2]
+ vertical_zmax_box=cubit.get_bounding_box('surface',k)[7]
+ normal=cubit.get_surface_normal(k)
+ check=(vertical_centerpoint >= vertical_centerpoint_top) and (vertical_zmax_box >= vertical_zmax_box_top) and (normal >= normal_top)
+ if check:
+ top=k
+ if top in surf_vertical:
+ surf_vertical.remove(top)
+ if top not in surf_or:
+ surf_or.append(top)
+ #if more than one surf is on the top, I get all the surfaces that are in touch with top surface but not the vertical surfaces
+ surftop=list(cubit.get_adjacent_surfaces("surface", top)) #top is included in the list
for s in surf_vertical:
try:
surftop.remove(s)
except:
pass
top=surftop
+ #check that all the surf are Horizontal or vertical
+ surf_all=surf_vertical+surf_or
+ if len(surf_all)!=len(lsurf):
+ print 'not all the surf are horizontal or vertical, check the normals'
+ print 'list of surfaces: ',surf_all
+ print 'list of vertical surface',surf_vertical
+ print 'list of horizontal surface',surf_or
+
bottom=[bottom]
return surf_or,surf_vertical,list_curve_or,list_curve_vertical,bottom,top
-
+def list2str(l):
+ if not isinstance(l,list): l=list(l)
+ return ' '.join(str(x) for x in l)
+
+def highlight(ent,l):
+ txt=list2str(l)
+ txt='highlight '+ent+' '+txt
+ cubit.cmd(txt)
Modified: seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/volumes.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/volumes.py 2013-09-26 00:49:25 UTC (rev 22849)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/geocubitlib/volumes.py 2013-09-26 00:57:45 UTC (rev 22850)
@@ -37,13 +37,20 @@
import start as start
print'volume'
cfg = start.start_cfg(filename=filename)
+ print cfg
#
if cfg.volume_type == 'layercake_volume_ascii_regulargrid_regularmap':
layercake_volume_ascii_regulargrid_mpiregularmap(filename=filename)
elif cfg.volume_type == 'layercake_volume_fromacis_mpiregularmap':
layercake_volume_fromacis_mpiregularmap(filename=filename)
+ elif cfg.volume_type == 'verticalsandwich_volume_ascii_regulargrid_mpiregularmap':
+ verticalsandwich_volume_ascii_regulargrid_mpiregularmap(filename=None)
+
-def layercake_volume_ascii_regulargrid_mpiregularmap(filename=None):
+
+
+
+def layercake_volume_ascii_regulargrid_mpiregularmap(filename=None,verticalsandwich=False):
import sys
import start as start
#
@@ -52,7 +59,7 @@
numpy = start.start_numpy()
cfg = start.start_cfg(filename=filename)
- from utilities import geo2utm, savegeometry,savesurf
+ from utilities import geo2utm, savegeometry,savesurf,cubit_command_check
from math import sqrt
#
@@ -63,21 +70,13 @@
#
#
command = "comment '"+"PROC: "+str(iproc)+"/"+str(numproc)+" '"
- cubit.cmd(command)
- #fer=open('ERROR_VOLUME_'+str(iproc),'w')
-
+ cubit_command_check(iproc,command,stop=True)
+ if verticalsandwich: cubit.cmd("comment 'Starting Vertical Sandwich'")
#
+ #get icpuy,icpux values
if mpiflag:
- x_slice=numpy.zeros([numproc],int)
- y_slice=numpy.zeros([numproc],int)
- for icpuy in range(0,cfg.nproc_eta):
- for icpux in range (0,cfg.nproc_xi):
- iprocnum=icpuy*cfg.nproc_xi+icpux
- #print iprocnum,cfg.nproc_xi,icpux,icpuy,cfg.nproc_xi,cfg.nproc_eta
- x_slice[iprocnum]=icpux
- y_slice[iprocnum]=icpuy
- icpux=x_slice[iproc]
- icpuy=y_slice[iproc]
+ icpux = iproc % cfg.nproc_xi
+ icpuy = int(iproc / cfg.nproc_xi)
else:
icpuy=int(cfg.id_proc/cfg.nproc_xi)
icpux=cfg.id_proc%cfg.nproc_xi
@@ -157,7 +156,7 @@
#print coordx,coordy,nx,ny
#
print 'end of building grid '+str(iproc)
- print 'number of point: ', len(coordx)
+ print 'number of point: ', len(coordx)*len(coordy)
#
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#for each processor
@@ -180,6 +179,17 @@
#
#create vertex
for inz in range(0,cfg.nz):
+ if cfg.sea and inz==cfg.nz-1: #sea layer
+ sealevel=True
+ bathymetry=False
+ elif cfg.sea and inz==cfg.nz-2: #bathymetry layer
+ sealevel=False
+ bathymetry=True
+ else:
+ sealevel=False
+ bathymetry=False
+ print sealevel,bathymetry
+
if cfg.bottomflat and inz == 0: #bottom layer
#
if cfg.geometry_format == 'ascii':
@@ -262,6 +272,7 @@
ivx=0
for ix in range(nxmin_cpu,nxmax_cpu+1):
zvertex=elev[ix,iy,inz]
+ #zvertex=adjust_sea_layers(zvertex,sealevel,bathymetry,cfg)
x_current,y_current=(coordx[ix,iy],coordy[ix,iy])
#
vertex.append(' Position '+ str( x_current ) +' '+ str( y_current )+' '+ str( zvertex ) )
@@ -323,32 +334,37 @@
#
#
#!create volume
- if cfg.osystem == 'macosx':
- pass
- elif cfg.osystem == 'linux':
- for inz in range(1,cfg.nz):
+ if not cfg.debugsurface:
+ if cfg.osystem == 'macosx':
+ pass
+ elif cfg.osystem == 'linux':
+ if cfg.nz == 1:
+ nsurface=2
+ else:
+ nsurface=cfg.nz
+ for inz in range(1,nsurface):
+ ner=cubit.get_error_count()
+ cubitcommand= 'create volume loft surface '+ str( inz+1 )+' '+str( inz )
+ cubit.cmd(cubitcommand)
+ ner2=cubit.get_error_count()
+ isurf=isurf+6
+ if ner == ner2:
+ cubitcommand= 'del surface 1 to '+ str( cfg.nz )
+ cubit.cmd(cubitcommand)
+ list_vol=cubit.parse_cubit_list("volume","all")
+ if len(list_vol) > 1:
+ cubitcommand= 'imprint volume all'
+ cubit.cmd(cubitcommand)
+ cubitcommand= 'merge all'
+ cubit.cmd(cubitcommand)
ner=cubit.get_error_count()
- cubitcommand= 'create volume loft surface '+ str( inz+1 )+' '+str( inz )
- cubit.cmd(cubitcommand)
- ner2=cubit.get_error_count()
- isurf=isurf+6
- if ner == ner2:
- cubitcommand= 'del surface 1 to '+ str( cfg.nz )
- cubit.cmd(cubitcommand)
- list_vol=cubit.parse_cubit_list("volume","all")
- if len(list_vol) > 1:
- cubitcommand= 'imprint volume all'
- cubit.cmd(cubitcommand)
- cubitcommand= 'merge all'
- cubit.cmd(cubitcommand)
- ner=cubit.get_error_count()
- #cubitcommand= 'composite create curve in vol all'
- #cubit.cmd(cubitcommand)
+ #cubitcommand= 'composite create curve in vol all'
+ #cubit.cmd(cubitcommand)
savegeometry(iproc,filename=filename)
- if cfg.geological_imprint:
- curvesname=[cfg.outlinebasin_curve,cfg.transition_curve,cfg.faulttrace_curve]
- outdir=cfg.working_dir
- imprint_topography_with_geological_outline(curvesname,outdir)
+ #if cfg.geological_imprint:
+ # curvesname=[cfg.outlinebasin_curve,cfg.transition_curve,cfg.faulttrace_curve]
+ # outdir=cfg.working_dir
+ # imprint_topography_with_geological_outline(curvesname,outdir)
#
#
cubit.cmd("set info "+cfg.cubit_info)
@@ -362,7 +378,6 @@
#
mpiflag,iproc,numproc,mpi = start.start_mpi()
#
- numpy = start.start_numpy()
cfg = start.start_cfg(filename=filename)
#
from utilities import geo2utm, savegeometry
@@ -409,16 +424,8 @@
cubit.cmd('move surface all x '+str(xmin)+' y '+str(ymin))
if mpiflag:
- x_slice=numpy.zeros([numproc],int)
- y_slice=numpy.zeros([numproc],int)
- for icpuy in range(0,cfg.nproc_eta):
- for icpux in range (0,cfg.nproc_xi):
- iprocnum=icpuy*cfg.nproc_xi+icpux
- #print iprocnum,cfg.nproc_xi,icpux,icpuy,cfg.nproc_xi,cfg.nproc_eta
- x_slice[iprocnum]=icpux
- y_slice[iprocnum]=icpuy
- icpux=x_slice[iproc]
- icpuy=y_slice[iproc]
+ icpux = iproc % cfg.nproc_xi
+ icpuy = int(iproc / cfg.nproc_xi)
else:
icpuy=int(cfg.id_proc/cfg.nproc_xi)
icpux=cfg.id_proc%cfg.nproc_xi
@@ -508,44 +515,45 @@
cubit.cmd(cubitcommand)
#
savegeometry(iproc,filename=filename)
- if cfg.geological_imprint:
- curvesname=[cfg.outlinebasin_curve,cfg.transition_curve,cfg.faulttrace_curve]
- outdir=cfg.working_dir
- imprint_topography_with_geological_outline(curvesname,outdir)
+ #if cfg.geological_imprint:
+ # curvesname=[cfg.outlinebasin_curve,cfg.transition_curve,cfg.faulttrace_curve]
+ # outdir=cfg.working_dir
+ # imprint_topography_with_geological_outline(curvesname,outdir)
-def imprint_topography_with_geological_outline(curvesname,outdir='.'):
- import sys,os
- from sets import Set
- #
- from utilities import load_curves,project_curves,get_v_h_list
-
- list_vol=cubit.parse_cubit_list("volume","all")
- surf_or,surf_vertical,list_curve_or,list_curve_vertical,bottom,top=get_v_h_list(list_vol)
-
-
-
- outlinebasin_curve=load_curves(curvesname[0])
- transition_curve=load_curves(curvesname[1])
- faulttrace_curve=load_curves(curvesname[2])
-
- curves=[]
- if outlinebasin_curve: curves=curves+outlinebasin_curve
- if transition_curve: curves=curves+transition_curve
- if faulttrace_curve: curves=curves+faulttrace_curve
-
- if curves:
- command='imprint tolerant surface '+str(top)+' with curve '+' '.join(str(x) for x in curves)+' merge'
- cubit.cmd(command)
-
-
- command = "merge surf all"
- cubit.cmd(command)
- command = "compress vol all"
- cubit.cmd(command)
- command = "compress surf all"
- cubit.cmd(command)
- command = "save as '"+outdirs+"/"+"imprinted_vol_"+str(iproc)+".cub' overwrite"
- cubit.cmd(command)
-
+#def imprint_topography_with_geological_outline(curvesname,outdir='.'):
+# import sys,os
+# from sets import Set
+# #
+# from utilities import load_curves,project_curves,get_v_h_list
+#
+# list_vol=cubit.parse_cubit_list("volume","all")
+# surf_or,surf_vertical,list_curve_or,list_curve_vertical,bottom,top=get_v_h_list(list_vol)
+#
+#
+#
+# outlinebasin_curve=load_curves(curvesname[0])
+# transition_curve=load_curves(curvesname[1])
+# faulttrace_curve=load_curves(curvesname[2])
+#
+# curves=[]
+# if outlinebasin_curve: curves=curves+outlinebasin_curve
+# if transition_curve: curves=curves+transition_curve
+# if faulttrace_curve: curves=curves+faulttrace_curve
+#
+# if curves:
+# command='imprint tolerant surface '+str(top)+' with curve '+' '.join(str(x) for x in curves)+' merge'
+# cubit.cmd(command)
+#
+#
+# command = "merge surf all"
+# cubit.cmd(command)
+# command = "compress vol all"
+# cubit.cmd(command)
+# command = "compress surf all"
+# cubit.cmd(command)
+# command = "save as '"+outdirs+"/"+"imprinted_vol_"+str(iproc)+".cub' overwrite"
+# cubit.cmd(command)
+#
+#
Deleted: seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/run_cubit2specfem3d.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/run_cubit2specfem3d.py 2013-09-26 00:49:25 UTC (rev 22849)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/run_cubit2specfem3d.py 2013-09-26 00:57:45 UTC (rev 22850)
@@ -1,40 +0,0 @@
-#!/usr/bin/env python
-
-import os
-import sys
-
-########################################################################################
-# put the path to your installation of the "bin" directory of the CUBIT package here
-########################################################################################
-sys.path.append("/home/komatits/bin/cubit/bin")
-import cubit
-
-# let us also add the local GEOCUBIT library
-sys.path.append("./geocubitlib")
-import boundary_definition
-import cubit2specfem3d
-
-########################################################################################
-# define the name of the CUBIT file to convert to SPECFEM format in the line below
-########################################################################################
-cubit.cmd('open "test_cpml_2layers.cub"')
-
-###### This is boundary_definition.py of GEOCUBIT
-#..... which extracts the bounding faces and defines them into blocks
-#reload(boundary_definition)
-
-# for CUBIT version 14 and higher (now called TRELIS) you may have to change 'face' to 'Face' in the line below
-boundary_definition.entities=['face']
-
-# this command can run in parallel or not
-#boundary_definition.define_bc(boundary_definition.entities,parallel=True)
-boundary_definition.define_bc(boundary_definition.entities,parallel=False)
-
-#### Export to SPECFEM3D format using cubit2specfem3d.py of GEOCUBIT
-
-os.system('mkdir -p MESH')
-#reload(cubit2specfem3d)
-cubit2specfem3d.export2SPECFEM3D('MESH')
-
-# all files needed by SCOTCH are now in directory MESH...
-
Modified: seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/setup.py
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/setup.py 2013-09-26 00:49:25 UTC (rev 22849)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/setup.py 2013-09-26 00:57:45 UTC (rev 22850)
@@ -6,6 +6,7 @@
description='CUBIT plugin',
author='emanuele casarotti',
author_email='emanuele.casarotti at ingv.it',
+ license='GNU General Public License v3 or later',
url='',
packages=['','geocubitlib'],
scripts=['GEOCUBIT.py']
Deleted: seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/test_cpml_2layers.cub
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/test_cpml_2layers.cub 2013-09-26 00:49:25 UTC (rev 22849)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/test_cpml_2layers.cub 2013-09-26 00:57:45 UTC (rev 22850)
@@ -1 +0,0 @@
-link ../examples/CPML_examples/CPML_realistic_mesh_with_topography/test_cpml_2layers.cub
\ No newline at end of file
Deleted: seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/utility/jq.sh
===================================================================
--- seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/utility/jq.sh 2013-09-26 00:49:25 UTC (rev 22849)
+++ seismo/3D/SPECFEM3D/trunk/CUBIT_GEOCUBIT/utility/jq.sh 2013-09-26 00:57:45 UTC (rev 22850)
@@ -1,8 +0,0 @@
-#!/bin/sh -f
-#PBS -N GEOCUBIT
-#PBS -o joba.log
-#PBS -l walltime=2:0:0
-#PBS -m abe -M emanuele.casarotti at ingv.it
-#PBS -j oe
-cd /home/casarotti/geocubit4
-python GEOCUBIT.py --mesh --build_volume --cfg='/home/casarotti/geocubit4/example/abruzzo.cfg' --id_proc=${PBS_ARRAY_INDEX} > ${PBS_ARRAY_INDEX}.log
More information about the CIG-COMMITS
mailing list