diff --git a/.gitignore b/.gitignore deleted file mode 100644 index db7a2db..0000000 --- a/.gitignore +++ /dev/null @@ -1,40 +0,0 @@ -# Compiled Object files -*.slo -*.lo -*.o -*.obj - -# Precompiled Headers -*.gch -*.pch - -# Compiled Dynamic libraries -*.so -*.dylib -*.dll - -# Fortran module files -*.mod - -# Compiled Static libraries -*.lai -*.la -*.a -*.lib - -# Executables -*.exe -*.out -*.app - -#Backup files -*~ - -#Old files -*.old - -#Executables -bin/* - -#Build directory -build/* diff --git a/CMake/FindCBLAS.cmake b/CMake/FindCBLAS.cmake deleted file mode 100644 index b729bd7..0000000 --- a/CMake/FindCBLAS.cmake +++ /dev/null @@ -1,21 +0,0 @@ - -FIND_PATH(CBLAS_INCLUDE_DIR NAMES cblas.h HINTS ${CBLAS_INCLUDE_CUSTOM_PATHS} PATHS ${CBLAS_INCLUDE_SEARCH_PATHS}) -FIND_LIBRARY(CBLAS_LIBRARIES NAMES cblas HINTS ${CBLAS_LIB_CUSTOM_PATHS} PATHS ${CBLAS_LIB_SEARCH_PATHS}) - -FIND_LIBRARY(BLAS_LIBRARIES NAMES blas HINTS ${BLAS_LIB_CUSTOM_PATHS} PATHS ${BLAS_LIB_SEARCH_PATHS}) - -SET(LOOKED_FOR - CBLAS_INCLUDE_DIR - CBLAS_LIBRARIES - BLAS_LIBRARIES -) - -include(FindPackageHandleStandardArgs) -find_package_handle_standard_args(CBLAS DEFAULT_MSG ${LOOKED_FOR}) - -IF(CBLAS_FOUND) - MESSAGE(STATUS "Found CBLAS & BLAS") - MESSAGE(STATUS CBLAS_INCLUDE_DIR=${CBLAS_INCLUDE_DIR}) - MESSAGE(STATUS CBLAS_LIBRARIES=${CBLAS_LIBRARIES}) - MESSAGE(STATUS BLAS_LIBRARIES=${BLAS_LIBRARIES}) -ENDIF(CBLAS_FOUND) diff --git a/CMake/FindLAPACKE.cmake b/CMake/FindLAPACKE.cmake deleted file mode 100644 index 5b8bf5b..0000000 --- a/CMake/FindLAPACKE.cmake +++ /dev/null @@ -1,17 +0,0 @@ - -FIND_PATH(LAPACKE_INCLUDE_DIR NAMES lapacke.h HINTS ${LAPACKE_INCLUDE_CUSTOM_PATHS} PATHS ${LAPACKE_INCLUDE_SEARCH_PATHS}) -FIND_LIBRARY(LAPACKE_LIBRARIES NAMES lapacke HINTS ${LAPACKE_LIB_CUSTOM_PATHS} PATHS ${LAPACKE_LIB_SEARCH_PATHS}) - -SET(LOOKED_FOR - LAPACKE_INCLUDE_DIR - LAPACKE_LIBRARIES -) - -include(FindPackageHandleStandardArgs) -find_package_handle_standard_args(LAPACKE DEFAULT_MSG ${LOOKED_FOR}) - -IF(LAPACKE_FOUND) - MESSAGE(STATUS "Found LAPACKE") - MESSAGE(STATUS LAPACKE_INCLUDE_DIR=${LAPACKE_INCLUDE_DIR}) - MESSAGE(STATUS LAPACKE_LIBRARIES=${LAPACKE_LIBRARIES}) -ENDIF(LAPACKE_FOUND) diff --git a/CMake/ubuntu.cmake b/CMake/ubuntu.cmake deleted file mode 100644 index 33ebc43..0000000 --- a/CMake/ubuntu.cmake +++ /dev/null @@ -1,33 +0,0 @@ - -SET(CBLAS_INCLUDE_SEARCH_PATHS - /usr/include - /usr/include/atlas - /usr/include/openblas - CACHE PATH "" -) - -SET(CBLAS_LIB_SEARCH_PATHS - /usr/lib - /usr/lib/atlas-base - /usr/lib/openblas-base - CACHE PATH "" -) - -SET(BLAS_LIB_SEARCH_PATHS - /usr/lib - /usr/lib/libblas - /usr/lib/atlas-base - /usr/lib/atlas-base/atlas - /usr/lib/openblas-base - CACHE PATH "" -) - -SET(LAPACKE_INCLUDE_SEARCH_PATHS - /usr/include - CACHE PATH "" -) - -SET(LAPACKE_LIB_SEARCH_PATHS - /usr/lib - CACHE PATH "" -) diff --git a/CMakeLists.txt b/CMakeLists.txt deleted file mode 100644 index 87de427..0000000 --- a/CMakeLists.txt +++ /dev/null @@ -1,55 +0,0 @@ -# Main CMake file for NativeSolid project - -PROJECT(NativeSolid) -CMAKE_MINIMUM_REQUIRED(VERSION 2.6) - -SET(LIBRARY_OUTPUT_PATH ${PROJECT_BINARY_DIR}/../bin CACHE PATH - "Single output directory for building all libraries.") -SET(EXECUTABLE_OUTPUT_PATH ${PROJECT_BINARY_DIR}/../bin CACHE PATH - "Single output directory for building all executables.") - -set(CMAKE_MACOSX_RPATH ON) -set(CMAKE_SKIP_BUILD_RPATH FALSE) -set(CMAKE_BUILD_WITH_INSTALL_RPATH FALSE) - -MARK_AS_ADVANCED(LIBRARY_OUTPUT_PATH EXECUTABLE_OUTPUT_PATH) - -if( NOT CMAKE_BUILD_TYPE ) - set( CMAKE_BUILD_TYPE Release CACHE STRING - "Choose the type of build, options are: None Debug Release RelWithDebInfo -MinSizeRel." - FORCE ) -endif() - -LIST(APPEND CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/CMake") - -# -- C++11 -- -INCLUDE(CheckCXXCompilerFlag) -CHECK_CXX_COMPILER_FLAG("-std=c++11" COMPILER_SUPPORTS_CXX11) -CHECK_CXX_COMPILER_FLAG("-std=c++0x" COMPILER_SUPPORTS_CXX0X) -IF(COMPILER_SUPPORTS_CXX11) - SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11") -ELSEIF(COMPILER_SUPPORTS_CXX0X) - SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++0x") -ELSE() - MESSAGE(STATUS "The compiler ${CMAKE_CXX_COMPILER} has no C++11 support. Please use a different C++ compiler.") -ENDIF() - -# -- CBLAS/LAPACKE -- -FIND_PACKAGE(LAPACKE REQUIRED) -FIND_PACKAGE(CBLAS REQUIRED) - -# -- Python -- -FIND_PACKAGE(PythonLibs REQUIRED) - MESSAGE(STATUS "PYTHON_INCLUDE_PATH=${PYTHON_INCLUDE_PATH}") - -# -- SWIG -- -FIND_PACKAGE(SWIG REQUIRED) -SET(CMAKE_SWIG_OUTDIR "${EXECUTABLE_OUTPUT_PATH}") - -INCLUDE_DIRECTORIES(${PROJECT_SOURCE_DIR}/include ${PROJECT_SOURCE_DIR}/api ${LAPACKE_INCLUDE_DIR} ${CBLAS_INCLUDE_DIR} ${PYTHON_INCLUDE_PATH}) - -# -- Sub directories -- -ADD_SUBDIRECTORY( src ) -ADD_SUBDIRECTORY( _api ) - diff --git a/INSTALL.txt b/INSTALL.txt deleted file mode 100644 index b898470..0000000 --- a/INSTALL.txt +++ /dev/null @@ -1,24 +0,0 @@ -Installation of NativeSolid code - -Tools required : -*CMAKE (~sudo apt-get install cmake) -*LAPACK/LAPACKE (liblapacke.so & lapacke.h) -*CBLAS (libcblas.a & cblas.h) -*C/C++ compiler (GCC) -*MPI implementatin (OpenMPI) - -1. Go to the CMake folder and set the right paths for headers and libraries in the FindCBLAS and FindLAPACKE files -2. Create and go to the build folder - $ mkdir build - $ cd build - and check that it is empty. If not : - $ rm -r * -3. Execute : - $ cmake .. - for configuration. -4. If no error, execute : - $ make (you can do : $ make -j numthreads , for multi-threads compilation) -5. Set your environment variable in you .bashrc file as : - export NATIVE_RUN="INSTALL_PATH/NativeSolid/bin" - export PATH=$PATH:$NATIVE_RUN - export PYTHONPATH=$PYTHONPATH:$NATIVE_RUN diff --git a/LICENSE b/LICENSE deleted file mode 100644 index 5f2dd7f..0000000 --- a/LICENSE +++ /dev/null @@ -1,505 +0,0 @@ - GNU LESSER GENERAL PUBLIC LICENSE - Version 2.1, February 1999 - - Copyright (C) 1991, 1999 Free Software Foundation, Inc. - 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA - Everyone is permitted to copy and distribute verbatim copies - of this license document, but changing it is not allowed. - -(This is the first released version of the Lesser GPL. It also counts - as the successor of the GNU Library Public License, version 2, hence - the version number 2.1.) - - Preamble - - The licenses for most software are designed to take away your -freedom to share and change it. By contrast, the GNU General Public -Licenses are intended to guarantee your freedom to share and change -free software--to make sure the software is free for all its users. - - This license, the Lesser General Public License, applies to some -specially designated software packages--typically libraries--of the -Free Software Foundation and other authors who decide to use it. You -can use it too, but we suggest you first think carefully about whether -this license or the ordinary General Public License is the better -strategy to use in any particular case, based on the explanations below. - - When we speak of free software, we are referring to freedom of use, -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 this service if you wish); that you receive source code or can get -it if you want it; that you can change the software and use pieces of -it in new free programs; and that you are informed that you can do -these things. - - To protect your rights, we need to make restrictions that forbid -distributors to deny you these rights or to ask you to surrender these -rights. These restrictions translate to certain responsibilities for -you if you distribute copies of the library or if you modify it. - - For example, if you distribute copies of the library, whether gratis -or for a fee, you must give the recipients all the rights that we gave -you. You must make sure that they, too, receive or can get the source -code. If you link other code with the library, you must provide -complete object files to the recipients, so that they can relink them -with the library after making changes to the library and recompiling -it. And you must show them these terms so they know their rights. - - We protect your rights with a two-step method: (1) we copyright the -library, and (2) we offer you this license, which gives you legal -permission to copy, distribute and/or modify the library. - - To protect each distributor, we want to make it very clear that -there is no warranty for the free library. Also, if the library is -modified by someone else and passed on, the recipients should know -that what they have is not the original version, so that the original -author's reputation will not be affected by problems that might be -introduced by others. - - Finally, software patents pose a constant threat to the existence of -any free program. We wish to make sure that a company cannot -effectively restrict the users of a free program by obtaining a -restrictive license from a patent holder. Therefore, we insist that -any patent license obtained for a version of the library must be -consistent with the full freedom of use specified in this license. - - Most GNU software, including some libraries, is covered by the -ordinary GNU General Public License. This license, the GNU Lesser -General Public License, applies to certain designated libraries, and -is quite different from the ordinary General Public License. We use -this license for certain libraries in order to permit linking those -libraries into non-free programs. - - When a program is linked with a library, whether statically or using -a shared library, the combination of the two is legally speaking a -combined work, a derivative of the original library. The ordinary -General Public License therefore permits such linking only if the -entire combination fits its criteria of freedom. The Lesser General -Public License permits more lax criteria for linking other code with -the library. - - We call this license the "Lesser" General Public License because it -does Less to protect the user's freedom than the ordinary General -Public License. It also provides other free software developers Less -of an advantage over competing non-free programs. These disadvantages -are the reason we use the ordinary General Public License for many -libraries. However, the Lesser license provides advantages in certain -special circumstances. - - For example, on rare occasions, there may be a special need to -encourage the widest possible use of a certain library, so that it becomes -a de-facto standard. To achieve this, non-free programs must be -allowed to use the library. A more frequent case is that a free -library does the same job as widely used non-free libraries. In this -case, there is little to gain by limiting the free library to free -software only, so we use the Lesser General Public License. - - In other cases, permission to use a particular library in non-free -programs enables a greater number of people to use a large body of -free software. For example, permission to use the GNU C Library in -non-free programs enables many more people to use the whole GNU -operating system, as well as its variant, the GNU/Linux operating -system. - - Although the Lesser General Public License is Less protective of the -users' freedom, it does ensure that the user of a program that is -linked with the Library has the freedom and the wherewithal to run -that program using a modified version of the Library. - - The precise terms and conditions for copying, distribution and -modification follow. Pay close attention to the difference between a -"work based on the library" and a "work that uses the library". The -former contains code derived from the library, whereas the latter must -be combined with the library in order to run. - - GNU LESSER GENERAL PUBLIC LICENSE - TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION - - 0. This License Agreement applies to any software library or other -program which contains a notice placed by the copyright holder or -other authorized party saying it may be distributed under the terms of -this Lesser General Public License (also called "this License"). -Each licensee is addressed as "you". - - A "library" means a collection of software functions and/or data -prepared so as to be conveniently linked with application programs -(which use some of those functions and data) to form executables. - - The "Library", below, refers to any such software library or work -which has been distributed under these terms. A "work based on the -Library" means either the Library or any derivative work under -copyright law: that is to say, a work containing the Library or a -portion of it, either verbatim or with modifications and/or translated -straightforwardly into another language. (Hereinafter, translation is -included without limitation in the term "modification".) - - "Source code" for a work means the preferred form of the work for -making modifications to it. For a library, complete source code means -all the source code for all modules it contains, plus any associated -interface definition files, plus the scripts used to control compilation -and installation of the library. - - Activities other than copying, distribution and modification are not -covered by this License; they are outside its scope. The act of -running a program using the Library is not restricted, and output from -such a program is covered only if its contents constitute a work based -on the Library (independent of the use of the Library in a tool for -writing it). Whether that is true depends on what the Library does -and what the program that uses the Library does. - - 1. You may copy and distribute verbatim copies of the Library's -complete source code as you receive it, in any medium, provided that -you conspicuously and appropriately publish on each copy an -appropriate copyright notice and disclaimer of warranty; keep intact -all the notices that refer to this License and to the absence of any -warranty; and distribute a copy of this License along with the -Library. - - You may charge a fee for the physical act of transferring a copy, -and you may at your option offer warranty protection in exchange for a -fee. - - 2. You may modify your copy or copies of the Library or any portion -of it, thus forming a work based on the Library, and copy and -distribute such modifications or work under the terms of Section 1 -above, provided that you also meet all of these conditions: - - a) The modified work must itself be a software library. - - b) You must cause the files modified to carry prominent notices - stating that you changed the files and the date of any change. - - c) You must cause the whole of the work to be licensed at no - charge to all third parties under the terms of this License. - - d) If a facility in the modified Library refers to a function or a - table of data to be supplied by an application program that uses - the facility, other than as an argument passed when the facility - is invoked, then you must make a good faith effort to ensure that, - in the event an application does not supply such function or - table, the facility still operates, and performs whatever part of - its purpose remains meaningful. - - (For example, a function in a library to compute square roots has - a purpose that is entirely well-defined independent of the - application. Therefore, Subsection 2d requires that any - application-supplied function or table used by this function must - be optional: if the application does not supply it, the square - root function must still compute square roots.) - -These requirements apply to the modified work as a whole. If -identifiable sections of that work are not derived from the Library, -and can be reasonably considered independent and separate works in -themselves, then this License, and its terms, do not apply to those -sections when you distribute them as separate works. But when you -distribute the same sections as part of a whole which is a work based -on the Library, the distribution of the whole must be on the terms of -this License, whose permissions for other licensees extend to the -entire whole, and thus to each and every part regardless of who wrote -it. - -Thus, it is not the intent of this section to claim rights or contest -your rights to work written entirely by you; rather, the intent is to -exercise the right to control the distribution of derivative or -collective works based on the Library. - -In addition, mere aggregation of another work not based on the Library -with the Library (or with a work based on the Library) on a volume of -a storage or distribution medium does not bring the other work under -the scope of this License. - - 3. You may opt to apply the terms of the ordinary GNU General Public -License instead of this License to a given copy of the Library. To do -this, you must alter all the notices that refer to this License, so -that they refer to the ordinary GNU General Public License, version 2, -instead of to this License. (If a newer version than version 2 of the -ordinary GNU General Public License has appeared, then you can specify -that version instead if you wish.) Do not make any other change in -these notices. - - Once this change is made in a given copy, it is irreversible for -that copy, so the ordinary GNU General Public License applies to all -subsequent copies and derivative works made from that copy. - - This option is useful when you wish to copy part of the code of -the Library into a program that is not a library. - - 4. You may copy and distribute the Library (or a portion or -derivative of it, under Section 2) in object code or executable form -under the terms of Sections 1 and 2 above provided that you accompany -it with the complete corresponding machine-readable source code, which -must be distributed under the terms of Sections 1 and 2 above on a -medium customarily used for software interchange. - - If distribution of object code is made by offering access to copy -from a designated place, then offering equivalent access to copy the -source code from the same place satisfies the requirement to -distribute the source code, even though third parties are not -compelled to copy the source along with the object code. - - 5. A program that contains no derivative of any portion of the -Library, but is designed to work with the Library by being compiled or -linked with it, is called a "work that uses the Library". Such a -work, in isolation, is not a derivative work of the Library, and -therefore falls outside the scope of this License. - - However, linking a "work that uses the Library" with the Library -creates an executable that is a derivative of the Library (because it -contains portions of the Library), rather than a "work that uses the -library". The executable is therefore covered by this License. -Section 6 states terms for distribution of such executables. - - When a "work that uses the Library" uses material from a header file -that is part of the Library, the object code for the work may be a -derivative work of the Library even though the source code is not. -Whether this is true is especially significant if the work can be -linked without the Library, or if the work is itself a library. The -threshold for this to be true is not precisely defined by law. - - If such an object file uses only numerical parameters, data -structure layouts and accessors, and small macros and small inline -functions (ten lines or less in length), then the use of the object -file is unrestricted, regardless of whether it is legally a derivative -work. (Executables containing this object code plus portions of the -Library will still fall under Section 6.) - - Otherwise, if the work is a derivative of the Library, you may -distribute the object code for the work under the terms of Section 6. -Any executables containing that work also fall under Section 6, -whether or not they are linked directly with the Library itself. - - 6. As an exception to the Sections above, you may also combine or -link a "work that uses the Library" with the Library to produce a -work containing portions of the Library, and distribute that work -under terms of your choice, provided that the terms permit -modification of the work for the customer's own use and reverse -engineering for debugging such modifications. - - You must give prominent notice with each copy of the work that the -Library is used in it and that the Library and its use are covered by -this License. You must supply a copy of this License. If the work -during execution displays copyright notices, you must include the -copyright notice for the Library among them, as well as a reference -directing the user to the copy of this License. Also, you must do one -of these things: - - a) Accompany the work with the complete corresponding - machine-readable source code for the Library including whatever - changes were used in the work (which must be distributed under - Sections 1 and 2 above); and, if the work is an executable linked - with the Library, with the complete machine-readable "work that - uses the Library", as object code and/or source code, so that the - user can modify the Library and then relink to produce a modified - executable containing the modified Library. (It is understood - that the user who changes the contents of definitions files in the - Library will not necessarily be able to recompile the application - to use the modified definitions.) - - b) Use a suitable shared library mechanism for linking with the - Library. A suitable mechanism is one that (1) uses at run time a - copy of the library already present on the user's computer system, - rather than copying library functions into the executable, and (2) - will operate properly with a modified version of the library, if - the user installs one, as long as the modified version is - interface-compatible with the version that the work was made with. - - c) Accompany the work with a written offer, valid for at - least three years, to give the same user the materials - specified in Subsection 6a, above, for a charge no more - than the cost of performing this distribution. - - d) If distribution of the work is made by offering access to copy - from a designated place, offer equivalent access to copy the above - specified materials from the same place. - - e) Verify that the user has already received a copy of these - materials or that you have already sent this user a copy. - - For an executable, the required form of the "work that uses the -Library" must include any data and utility programs needed for -reproducing the executable from it. However, as a special exception, -the materials to be distributed need not include anything that is -normally distributed (in either source or binary form) with the major -components (compiler, kernel, and so on) of the operating system on -which the executable runs, unless that component itself accompanies -the executable. - - It may happen that this requirement contradicts the license -restrictions of other proprietary libraries that do not normally -accompany the operating system. Such a contradiction means you cannot -use both them and the Library together in an executable that you -distribute. - - 7. You may place library facilities that are a work based on the -Library side-by-side in a single library together with other library -facilities not covered by this License, and distribute such a combined -library, provided that the separate distribution of the work based on -the Library and of the other library facilities is otherwise -permitted, and provided that you do these two things: - - a) Accompany the combined library with a copy of the same work - based on the Library, uncombined with any other library - facilities. This must be distributed under the terms of the - Sections above. - - b) Give prominent notice with the combined library of the fact - that part of it is a work based on the Library, and explaining - where to find the accompanying uncombined form of the same work. - - 8. You may not copy, modify, sublicense, link with, or distribute -the Library except as expressly provided under this License. Any -attempt otherwise to copy, modify, sublicense, link with, or -distribute the Library is void, and will automatically terminate your -rights under this License. However, parties who have received copies, -or rights, from you under this License will not have their licenses -terminated so long as such parties remain in full compliance. - - 9. You are not required to accept this License, since you have not -signed it. However, nothing else grants you permission to modify or -distribute the Library or its derivative works. These actions are -prohibited by law if you do not accept this License. Therefore, by -modifying or distributing the Library (or any work based on the -Library), you indicate your acceptance of this License to do so, and -all its terms and conditions for copying, distributing or modifying -the Library or works based on it. - - 10. Each time you redistribute the Library (or any work based on the -Library), the recipient automatically receives a license from the -original licensor to copy, distribute, link with or modify the Library -subject to these terms and conditions. You may not impose any further -restrictions on the recipients' exercise of the rights granted herein. -You are not responsible for enforcing compliance by third parties with -this License. - - 11. If, as a consequence of a court judgment or allegation of patent -infringement or for any other reason (not limited to patent issues), -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 -distribute so as to satisfy simultaneously your obligations under this -License and any other pertinent obligations, then as a consequence you -may not distribute the Library at all. For example, if a patent -license would not permit royalty-free redistribution of the Library by -all those who receive copies directly or indirectly through you, then -the only way you could satisfy both it and this License would be to -refrain entirely from distribution of the Library. - -If any portion of this section is held invalid or unenforceable under any -particular circumstance, the balance of the section is intended to apply, -and the section as a whole is intended to apply in other circumstances. - -It is not the purpose of this section to induce you to infringe any -patents or other property right claims or to contest validity of any -such claims; this section has the sole purpose of protecting the -integrity of the free software distribution system which is -implemented by public license practices. Many people have made -generous contributions to the wide range of software distributed -through that system in reliance on consistent application of that -system; it is up to the author/donor to decide if he or she is willing -to distribute software through any other system and a licensee cannot -impose that choice. - -This section is intended to make thoroughly clear what is believed to -be a consequence of the rest of this License. - - 12. If the distribution and/or use of the Library is restricted in -certain countries either by patents or by copyrighted interfaces, the -original copyright holder who places the Library under this License may add -an explicit geographical distribution limitation excluding those countries, -so that distribution is permitted only in or among countries not thus -excluded. In such case, this License incorporates the limitation as if -written in the body of this License. - - 13. The Free Software Foundation may publish revised and/or new -versions of the Lesser 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 Library -specifies a version number of this License which applies to it and -"any later version", you have the option of following the terms and -conditions either of that version or of any later version published by -the Free Software Foundation. If the Library does not specify a -license version number, you may choose any version ever published by -the Free Software Foundation. - - 14. If you wish to incorporate parts of the Library into other free -programs whose distribution conditions are incompatible with these, -write to the author to ask for permission. For software which is -copyrighted by the Free Software Foundation, write to the Free -Software Foundation; we sometimes make exceptions for this. Our -decision will be guided by the two goals of preserving the free status -of all derivatives of our free software and of promoting the sharing -and reuse of software generally. - - NO WARRANTY - - 15. BECAUSE THE LIBRARY IS LICENSED FREE OF CHARGE, THERE IS NO -WARRANTY FOR THE LIBRARY, TO THE EXTENT PERMITTED BY APPLICABLE LAW. -EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR -OTHER PARTIES PROVIDE THE LIBRARY "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 -LIBRARY IS WITH YOU. SHOULD THE LIBRARY PROVE DEFECTIVE, YOU ASSUME -THE COST OF ALL NECESSARY SERVICING, REPAIR OR CORRECTION. - - 16. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN -WRITING WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY -AND/OR REDISTRIBUTE THE LIBRARY 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 -LIBRARY (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 LIBRARY TO OPERATE WITH ANY OTHER SOFTWARE), EVEN IF -SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH -DAMAGES. - - END OF TERMS AND CONDITIONS - - How to Apply These Terms to Your New Libraries - - If you develop a new library, and you want it to be of the greatest -possible use to the public, we recommend making it free software that -everyone can redistribute and change. You can do so by permitting -redistribution under these terms (or, alternatively, under the terms of the -ordinary General Public License). - - To apply these terms, attach the following notices to the library. It is -safest to attach them to the start of each source file to most effectively -convey the exclusion of warranty; and each file should have at least the -"copyright" line and a pointer to where the full notice is found. - - {description} - Copyright (C) {year} {fullname} - - This library is free software; you can redistribute it and/or - modify it under the terms of the GNU Lesser General Public - License as published by the Free Software Foundation; either - version 2.1 of the License, or (at your option) any later version. - - This library 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 - Lesser General Public License for more details. - - You should have received a copy of the GNU Lesser General Public - License along with this library; if not, write to the Free Software - Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 - USA - -Also add information on how to contact you by electronic and paper mail. - -You should also get your employer (if you work as a programmer) or your -school, if any, to sign a "copyright disclaimer" for the library, if -necessary. Here is a sample; alter the names: - - Yoyodyne, Inc., hereby disclaims all copyright interest in the - library `Frob' (a library for tweaking knobs) written by James Random - Hacker. - - {signature of Ty Coon}, 1 April 1990 - Ty Coon, President of Vice - -That's all there is to it! - diff --git a/RBMConf.cfg b/RBMConf.cfg deleted file mode 100644 index 4a12da3..0000000 --- a/RBMConf.cfg +++ /dev/null @@ -1,58 +0,0 @@ -% Unsteady (dynamic) simulation -UNSTEADY_SIMULATION = YES -% CSD solver (NATIVE or METAFOR) -CSD_SOLVER = NATIVE -% -MESH_FILE = 2D_FlatPlate_Rounded.su2 -% -MOVING_MARKER = plate -% -% If NATIVE -% Integration time step -DELTA_T = 0.001 -% Type of structural problem (SPRING_HOR or SPRING_VER or AIRFOIL) -STRUCT_TYPE = AIRFOIL -% Linearized equations of motion -LINEARIZE = YES -% Body mass -SPRING_MASS = 1.356 -%Inertia around CG -INERTIA_CG = 1000 -%Inertia around flexural axis -INERTIA_FLEXURAL = 2.07578e-4 -% Spring stiffness -SPRING_STIFFNESS = 2658.623 -% Spring damping -SPRING_DAMPING = 0.240172 -% -TORSIONAL_STIFFNESS = 0.663817 -% -TORSIONAL_DAMPING = 3.5216e-4 -% -CORD = 0.035 -%Position of the flexural axis -FLEXURAL_AXIS = -0.0028 -%Center of gravity -GRAVITY_CENTER = 0.0 -% -% INITIAL CONDITIONS -%Plunging -INITIAL_DISP = 0.0 -%Pitching -INITIAL_ANGLE = 0.7854 -% -% Restart solution from previous computations -RESTART_SOL = NO -% -DELTA_ITER_WRITE = 10 -% Restart solution file name -RESTART_FILE = restart_solid.dat -% Start time (usually 0) in seconds -START_TIME = 0 -% End time in seconds -STOP_TIME = 1.0 -% Integration algorithm (ALPHAGEN, RK4) -INTEGRATION_ALGO = ALPHAGEN -%INTEGRATION_ALGO = RK4 -% Integration parameter -RHO = 0.0 diff --git a/README.md b/README.md index 85359a6..9584408 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,3 @@ -# NativeSolid -Native solid solver for FSI +# NativeSolid / RBM + +This repository has been moved to the [GitLab instance of the University of Liège](https://gitlab.uliege.be/am-dept/NativeSolid) diff --git a/_api/CMakeLists.txt b/_api/CMakeLists.txt deleted file mode 100644 index f9a5495..0000000 --- a/_api/CMakeLists.txt +++ /dev/null @@ -1,19 +0,0 @@ -# CMake input file of the SWIG wrapper around "Native.so" - -INCLUDE(${SWIG_USE_FILE}) - -FILE(GLOB SRCS *.h *.cpp *.inl *.swg) -FILE(GLOB ISRCS *.i) - -SET(CMAKE_SWIG_FLAGS "") -SET_SOURCE_FILES_PROPERTIES(${ISRCS} PROPERTIES CPLUSPLUS ON) - -SET(SWINCFLAGS --I${PROJECT_SOURCE_DIR}/api -) - -SET_SOURCE_FILES_PROPERTIES(${ISRCS} PROPERTIES SWIG_FLAGS "${SWINCFLAGS}") - -SWIG_ADD_MODULE(NativeSolid python ${ISRCS} ${SRCS}) - -SWIG_LINK_LIBRARIES(NativeSolid Native ${PYTHON_LIBRARIES}) diff --git a/_api/Native.i b/_api/Native.i deleted file mode 100644 index 2ea59d2..0000000 --- a/_api/Native.i +++ /dev/null @@ -1,22 +0,0 @@ -// SWIG input file of the 'NativeSolid' module - -%feature("autodoc","1"); - -%module(docstring= -"'NativeSolid' module", -directors="1", -threads="1" -) NativeSolid -%{ - -#include -#include "NativeSolid_API.h" - - -%} - -// ----------- MODULES UTILISES ------------ -%include "std_string.i" - -// ----------- NATIVE CLASSES ---------------- -%include "NativeSolid_API.h" diff --git a/api/NativeSolid_API.cpp b/api/NativeSolid_API.cpp deleted file mode 100644 index 74e965f..0000000 --- a/api/NativeSolid_API.cpp +++ /dev/null @@ -1,893 +0,0 @@ -#include "NativeSolid_API.h" -#include - -using namespace std; - -const int MASTER_NODE = 0; - - -NativeSolidSolver::NativeSolidSolver(string str, bool FSIComp):confFile(str){ - - int rank = MASTER_NODE; - - historyFile.open("NativeHistory.dat", ios::out); - -/*--- MPI initialization, and buffer setting ---*/ -#ifdef HAVE_MPI - MPI_Comm_rank(MPI_COMM_WORLD, &rank); -#endif - - if(rank == MASTER_NODE){ - if(FSIComp)cout << endl << "***************************** Setting NativeSolid for FSI simulation *****************************" << endl; - else cout << endl <<"***************************** Setting NativeSolid for CSD simulation *****************************" << endl; - } - - config = NULL; - output = NULL; - geometry = NULL; - structure = NULL; - integrator = NULL; - - /*--- Initialize the main containers ---*/ - config = new Config(confFile); - output = new Output(); - - /*--- Read CSD configuration file ---*/ - cout << endl << "\n----------------------- Reading Native configuration file ----------------------" << endl; - config->ReadConfig(); - - /*--- Read a SU2 native mesh file ---*/ - cout << endl << "\n----------------------- Reading SU2 based mesh file ----------------------" << endl; - geometry = new Geometry(config); - - /*--- Initialize structural container and create the structural model ---*/ - cout << endl << "\n----------------------- Creating the structural model ----------------------" << endl; - structure = new Structure(config); - - cout << endl << "\n----------------------- Creating the FSI interface ----------------------" << endl; - double* Coord; - unsigned long iMarker(0), iPoint; - - while(iMarker < geometry->GetnMarkers()){ - if(geometry->GetMarkersMoving(iMarker)){ - nSolidInterfaceVertex = geometry->nVertex[iMarker]; - break; - } - iMarker++; - } - - varCoordNorm = 0.0; - - cout << nSolidInterfaceVertex << " nodes on the moving interface have to be tracked." << endl; - - /*--- Initialize the temporal integrator ---*/ - cout << endl << "\n----------------------- Setting integration parameter ----------------------" << endl; - integrator = new Integration(config, structure); - if(config->GetUnsteady() == "YES"){ - integrator->SetInitialConditions(config, structure); - } - - cout << endl << "\n----------------------- Setting FSI features ----------------------" << endl; - q_uM1.Initialize(structure->GetnDof()); - q_uM1.Reset(); - - if(rank == MASTER_NODE){ - if(structure->GetnDof() == 1){ - if(config->GetUnsteady() == "YES"){ - historyFile << "\"Time\"" << "\t" << "\"Displacement\"" << "\t" << "\"Velocity\"" << "\t" << "\"Acceleration\"" << endl; - } - else{ - historyFile << "\"FSI Iteration\"" << "\t" << "\"Displacement\"" << endl; - } - } - else if(structure->GetnDof() == 2){ - if(config->GetUnsteady() == "YES"){ - historyFile << "\"Time\"" << "\t" << "\"Displacement 1\"" << "\t" << "\"Displacement 2\"" << "\t" << "\"Velocity 1\"" << "\t" << "\"Velocity 2\"" << "\t" << "\"Acceleration 1\"" << "\t" << "\"Acceleration 2\"" << endl; - } - else{ - historyFile << "\"FSI Iteration\"" << "\t" << "\"Displacement 1\"" << "\t" << "\"Displacement 1\"" << endl; - } - } - } - - if(rank == MASTER_NODE){ - if(FSIComp)cout << endl << "***************************** NativeSolid is set for FSI simulation *****************************" << endl; - else cout << endl <<"***************************** NativeSolid is set for CSD simulation *****************************" << endl; - } - -} - -NativeSolidSolver::~NativeSolidSolver(){} - -void NativeSolidSolver::exit(){ - - int rank = MASTER_NODE; - -#ifdef HAVE_MPI - MPI_Comm_rank(MPI_COMM_WORLD, &rank); -#endif - - historyFile.close(); - if(rank == MASTER_NODE){ - cout << "Solid history file is closed." << endl; - cout << endl << "***************************** Exit NativeSolid *****************************" << endl; - } - - if (config != NULL) delete config; - if (geometry != NULL) delete geometry; - if (structure != NULL) delete structure; - if (integrator != NULL) delete integrator; - if (output != NULL) delete output; - -} - -void NativeSolidSolver::preprocessIteration(unsigned long ExtIter){ - - integrator->SetExtIter(ExtIter); -} - -void NativeSolidSolver::timeIteration(double t0, double tf){ - - integrator->TemporalIteration(t0, tf, structure); - - //mapRigidBodyMotion(false, false); - computeInterfacePosVel(false); - -} - -/* -void NativeSolidSolver::mapRigidBodyMotion(bool prediction, bool initialize){ - - double *Coord, *Coord_n, newCoord[3], Center[3], Center_n[3], newCenter[3], rotCoord[3], r[3]; - double varCoord[3] = {0.0, 0.0, 0.0}; - double rotMatrix[3][3] = {{0.0,0.0,0.0}, {0.0,0.0,0.0}, {0.0,0.0,0.0}}; - double dTheta, dPhi, dPsi; - double cosTheta, sinTheta, cosPhi, sinPhi, cosPsi, sinPsi; - unsigned short iMarker, iVertex, nDim(3); - unsigned long iPoint; - - double deltaT, q_n, qdot_n, qdot_nM1, q_nP1; - double alpha_n, alphadot_n, alphadot_nM1, alpha_nP1; - double alpha0, alpha1; - double disp, dAlpha; - double varCoordNorm2(0.0); - - //--- Get the current center of rotation (can vary at each iteration) --- - Center[0] = structure->GetCenterOfRotation_x(); - Center[1] = structure->GetCenterOfRotation_y(); - Center[2] = structure->GetCenterOfRotation_z(); - - //--- Get the center of rotation from previous time step --- - Center_n[0] = structure->GetCenterOfRotation_n_x(); - Center_n[1] = structure->GetCenterOfRotation_n_y(); - Center_n[2] = structure->GetCenterOfRotation_n_z(); - - if(!prediction){ - dTheta = 0.0; - dPhi = 0.0; - if (config->GetStructType() == "AIRFOIL"){ - dPsi = -( (*(integrator->GetDisp()))[1] - (*(integrator->GetDisp_n()))[1]); - newCenter[0] = Center[0]; - newCenter[1] = -(*(integrator->GetDisp()))[0]; - newCenter[2] = Center[2]; - disp = newCenter[1] - Center_n[1]; - } - else if(config->GetStructType() == "SPRING_HOR"){ - dPsi = 0.0; - newCenter[0] = (*(integrator->GetDisp()))[0]; - newCenter[1] = Center[1]; - newCenter[2] = Center[2]; - } - } - else{ - deltaT = config->GetDeltaT(); - alpha0 = 1.0; - alpha1 = 0.5; //Second order prediction - - q_n = (*(integrator->GetDisp()))[0]; - qdot_n = (*(integrator->GetVel()))[0]; - qdot_nM1 = (*(integrator->GetVel_n()))[0]; - q_nP1 = q_n + alpha0*deltaT*qdot_n + alpha1*deltaT*(qdot_n - qdot_nM1); - disp = q_nP1-q_n; - - if(structure->GetnDof() == 2){ - alpha_n = (*(integrator->GetDisp()))[1]; - alphadot_n = (*(integrator->GetVel()))[1]; - alphadot_nM1 = (*(integrator->GetVel_n()))[1]; - alpha_nP1 = alpha_n + alpha0*deltaT*alphadot_n + alpha1*deltaT*(alphadot_n - alphadot_nM1); - dAlpha = alpha_nP1-alpha_n; - } - dTheta = 0.0; - dPhi = 0.0; - if (config->GetStructType() == "AIRFOIL"){ - dPsi = -dAlpha; - newCenter[0] = Center[0]; - newCenter[1] = -q_nP1; - newCenter[2] = Center[2]; - } - else{ - dPsi = 0.0; - } - } - - cosTheta = cos(dTheta); cosPhi = cos(dPhi); cosPsi = cos(dPsi); - sinTheta = sin(dTheta); sinPhi = sin(dPhi); sinPsi = sin(dPsi); - - //--- Compute the rotation matrix. The implicit ordering is rotation about the x-axis, y-axis, then z-axis. --- - - rotMatrix[0][0] = cosPhi*cosPsi; - rotMatrix[1][0] = cosPhi*sinPsi; - rotMatrix[2][0] = -sinPhi; - - rotMatrix[0][1] = sinTheta*sinPhi*cosPsi - cosTheta*sinPsi; - rotMatrix[1][1] = sinTheta*sinPhi*sinPsi + cosTheta*cosPsi; - rotMatrix[2][1] = sinTheta*cosPhi; - - rotMatrix[0][2] = cosTheta*sinPhi*cosPsi + sinTheta*sinPsi; - rotMatrix[1][2] = cosTheta*sinPhi*sinPsi - sinTheta*cosPsi; - rotMatrix[2][2] = cosTheta*cosPhi; - - for(iMarker = 0; iMarker < geometry->GetnMarkers(); iMarker++){ - if (geometry->markersMoving[iMarker] == true){ - for(iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++){ - iPoint = geometry->vertex[iMarker][iVertex]; - Coord = geometry->node[iPoint]->GetCoord(); - Coord_n = geometry->node[iPoint]->GetCoord_n(); - - if(config->GetUnsteady() == "YES"){ - for (int iDim=0; iDim < nDim; iDim++){ - r[iDim] = Coord_n[iDim] - Center_n[iDim]; - } - } - else{ - for (int iDim=0; iDim < nDim; iDim++){ - r[iDim] = Coord[iDim] - Center[iDim]; - } - } - - rotCoord[0] = rotMatrix[0][0]*r[0] - + rotMatrix[0][1]*r[1] - + rotMatrix[0][2]*r[2]; - - rotCoord[1] = rotMatrix[1][0]*r[0] - + rotMatrix[1][1]*r[1] - + rotMatrix[1][2]*r[2]; - - rotCoord[2] = rotMatrix[2][0]*r[0] - + rotMatrix[2][1]*r[1] - + rotMatrix[2][2]*r[2]; - - for(int iDim=0; iDim < nDim; iDim++){ - newCoord[iDim] = newCenter[iDim] + rotCoord[iDim]; - varCoord[iDim] = newCoord[iDim] - Coord[iDim]; - } - - varCoordNorm2 += varCoord[0]*varCoord[0] + varCoord[1]*varCoord[1] + varCoord[2]*varCoord[2]; - - //--- Apply change of coordinates to the node on the moving interface --- - geometry->node[iPoint]->SetCoord(newCoord); - - //--- At initialisation, propagate the initial position of the inteface in the past --- - if(initialize) geometry->node[iPoint]->SetCoord_n(newCoord); - } - } - } - - varCoordNorm = sqrt(varCoordNorm2); - - //--- Update the position of the center of rotation --- - structure->SetCenterOfRotation_X(newCenter[0]); - structure->SetCenterOfRotation_Y(newCenter[1]); - structure->SetCenterOfRotation_Z(newCenter[2]); -} -*/ - -void NativeSolidSolver::computeInterfacePosVel(bool initialize){ - - double *Coord, *Coord_n, newCoord[3], newVel[3], Center[3], Center_n[3], newCenter[3], centerVel[3], rotCoord[3], r[3]; - double varCoord[3] = {0.0, 0.0, 0.0}; - double rotMatrix[3][3] = {{0.0,0.0,0.0}, {0.0,0.0,0.0}, {0.0,0.0,0.0}}; - double dTheta, dPhi, dPsi; - double psidot; - double cosTheta, sinTheta, cosPhi, sinPhi, cosPsi, sinPsi; - unsigned short iMarker, iVertex, nDim(3); - unsigned long iPoint; - double varCoordNorm2(0.0); - - /*--- Get the current center of rotation (can vary at each iteration) ---*/ - Center[0] = structure->GetCenterOfRotation_x(); - Center[1] = structure->GetCenterOfRotation_y(); - Center[2] = structure->GetCenterOfRotation_z(); - - /*--- Get the center of rotation from previous time step ---*/ - Center_n[0] = structure->GetCenterOfRotation_n_x(); - Center_n[1] = structure->GetCenterOfRotation_n_y(); - Center_n[2] = structure->GetCenterOfRotation_n_z(); - - dTheta = 0.0; - dPhi = 0.0; - if (config->GetStructType() == "AIRFOIL"){ - dPsi = -( (integrator->GetSolver()->GetDisp())[1] - (integrator->GetSolver()->GetDisp_n())[1] ); - psidot = (integrator->GetSolver()->GetVel())[1]; - newCenter[0] = Center[0]; - newCenter[1] = -(integrator->GetSolver()->GetDisp())[0]; - newCenter[2] = Center[2]; - centerVel[0] = 0.0; - centerVel[1] = -(integrator->GetSolver()->GetVel())[0]; - centerVel[2] = 0.0; - } - else if(config->GetStructType() == "SPRING_HOR"){ - dPsi = 0.0; - psidot = 0.0; - newCenter[0] = (integrator->GetSolver()->GetDisp())[0]; - newCenter[1] = Center[1]; - newCenter[2] = Center[2]; - centerVel[0] = (integrator->GetSolver()->GetVel())[0]; - centerVel[1] = 0.0; - centerVel[2] = 0.0; - } - else if(config->GetStructType() == "SPRING_VER"){ - dPsi = 0.0; - psidot = 0.0; - newCenter[0] = Center[0]; - newCenter[1] = (integrator->GetSolver()->GetDisp())[0]; - newCenter[2] = Center[2]; - centerVel[0] = 0.0; - centerVel[1] = (integrator->GetSolver()->GetVel())[0]; - centerVel[2] = 0.0; - } - else { } - - cosTheta = cos(dTheta); cosPhi = cos(dPhi); cosPsi = cos(dPsi); - sinTheta = sin(dTheta); sinPhi = sin(dPhi); sinPsi = sin(dPsi); - - /*--- Compute the rotation matrix. The implicit - ordering is rotation about the x-axis, y-axis, then z-axis. ---*/ - - rotMatrix[0][0] = cosPhi*cosPsi; - rotMatrix[1][0] = cosPhi*sinPsi; - rotMatrix[2][0] = -sinPhi; - - rotMatrix[0][1] = sinTheta*sinPhi*cosPsi - cosTheta*sinPsi; - rotMatrix[1][1] = sinTheta*sinPhi*sinPsi + cosTheta*cosPsi; - rotMatrix[2][1] = sinTheta*cosPhi; - - rotMatrix[0][2] = cosTheta*sinPhi*cosPsi + sinTheta*sinPsi; - rotMatrix[1][2] = cosTheta*sinPhi*sinPsi - sinTheta*cosPsi; - rotMatrix[2][2] = cosTheta*cosPhi; - - - for(iMarker = 0; iMarker < geometry->GetnMarkers(); iMarker++){ - if (geometry->markersMoving[iMarker] == true){ - for(iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++){ - - iPoint = geometry->vertex[iMarker][iVertex]; - Coord = geometry->node[iPoint]->GetCoord(); - Coord_n = geometry->node[iPoint]->GetCoord_n(); - - if(config->GetUnsteady() == "YES"){ - for (int iDim=0; iDim < nDim; iDim++){ - r[iDim] = Coord_n[iDim] - Center_n[iDim]; - } - } - else{ - for (int iDim=0; iDim < nDim; iDim++){ - r[iDim] = Coord[iDim] - Center[iDim]; - } - } - - rotCoord[0] = rotMatrix[0][0]*r[0] - + rotMatrix[0][1]*r[1] - + rotMatrix[0][2]*r[2]; - - rotCoord[1] = rotMatrix[1][0]*r[0] - + rotMatrix[1][1]*r[1] - + rotMatrix[1][2]*r[2]; - - rotCoord[2] = rotMatrix[2][0]*r[0] - + rotMatrix[2][1]*r[1] - + rotMatrix[2][2]*r[2]; - - for(int iDim=0; iDim < nDim; iDim++){ - newCoord[iDim] = newCenter[iDim] + rotCoord[iDim]; - varCoord[iDim] = newCoord[iDim] - Coord[iDim]; - } - - newVel[0] = centerVel[0] + psidot*(newCoord[1]-newCenter[1]); - newVel[1] = centerVel[1] - psidot*(newCoord[0]-newCenter[0]); - newVel[2] = centerVel[2] + 0.0; - - varCoordNorm2 += varCoord[0]*varCoord[0] + varCoord[1]*varCoord[1] + varCoord[2]*varCoord[2]; - - /*--- Apply change of coordinates to the node on the moving interface ---*/ - geometry->node[iPoint]->SetCoord(newCoord); - geometry->node[iPoint]->SetVel(newVel); - - /*--- At initialisation, propagate the initial position of the inteface in the past ---*/ - if(initialize){ - geometry->node[iPoint]->SetCoord_n(newCoord); - geometry->node[iPoint]->SetVel_n(newVel); - } - } - } - } - - varCoordNorm = sqrt(varCoordNorm2); - - /*--- Update the position of the center of rotation ---*/ - structure->SetCenterOfRotation_X(newCenter[0]); - structure->SetCenterOfRotation_Y(newCenter[1]); - structure->SetCenterOfRotation_Z(newCenter[2]); - -} - -void NativeSolidSolver::setInitialDisplacements(){ - //mapRigidBodyMotion(false, true); - computeInterfacePosVel(true); -} - -/* -void NativeSolidSolver::staticComputation(){ - - int rank = MASTER_NODE; - int size = 1; - -#ifdef HAVE_MPI - MPI_Comm_rank(MPI_COMM_WORLD, &rank); - MPI_Comm_size(MPI_COMM_WORLD, &size); -#endif - - integrator->StaticIteration(config,structure); - - if(rank == MASTER_NODE){ - if(structure->GetnDof() == 1){ - cout << "Static Displacement : " << (*(integrator->GetDisp()))[0] << endl; - } - else if(structure->GetnDof() == 2){ - cout << "Static Displacement 1 : " << (*(integrator->GetDisp()))[0] << endl; - cout << "Static Displacement 2 : " << (*(integrator->GetDisp()))[1] << endl; - } - } - - //mapRigidBodyMotion(false,false); - computeInterfacePosVel(false); - -} -*/ - -void NativeSolidSolver::writeSolution(double currentTime, double lastTime, double currentFSIIter, unsigned long ExtIter, unsigned long NbExtIter){ - - int rank = MASTER_NODE; - double DeltaT; - string restartFileName; - unsigned long DeltaIter = config->GetDeltaIterWrite(); - - DeltaT = currentTime-lastTime; - -#ifdef HAVE_MPI - MPI_Comm_rank(MPI_COMM_WORLD, &rank); -#endif - - if(rank == MASTER_NODE){ - if(structure->GetnDof() == 1){ - if(config->GetUnsteady() == "YES"){ - //if(currentTime == config->GetStartTime()) historyFile << "\"Time\"" << "\t" << "\"Displacement\"" << "\t" << "\"Velocity\"" << "\t" << "\"Acceleration\"" << "\t" << "\"Acceleration variable\"" << endl; - historyFile << currentTime << "\t" << (integrator->GetSolver()->GetDisp())[0] << "\t" << (integrator->GetSolver()->GetVel())[0] << "\t" << (integrator->GetSolver()->GetAcc())[0] << "\t" << (integrator->GetSolver()->GetAccVar())[0] << endl; - } - else{ - //if(currentFSIIter == config->GetStartTime()) historyFile << "\"FSI Iteration\"" << "\t" << "\"Displacement\"" << endl; - historyFile << currentFSIIter << "\t" << (integrator->GetSolver()->GetDisp())[0] << endl; - } - } - else if(structure->GetnDof() == 2){ - if(config->GetUnsteady() == "YES"){ - //if(currentTime == config->GetStartTime()) historyFile << "\"Time\"" << "\t" << "\"Displacement 1\"" << "\t" << "\"Displacement 2\"" << "\t" << "\"Velocity 1\"" << "\t" << "\"Velocity 2\"" << "\t" << "\"Acceleration 1\"" << "\t" << "\"Acceleration 2\"" << "\t" << "\"Acceleration variable 1\"" << "\t" << "\"Acceleration variable 2\"" << endl; - historyFile << currentTime << "\t" << (integrator->GetSolver()->GetDisp())[0] << "\t" << (integrator->GetSolver()->GetDisp())[1] << "\t" << (integrator->GetSolver()->GetVel())[0] << "\t" << (integrator->GetSolver()->GetVel())[1] << "\t" << (integrator->GetSolver()->GetAcc())[0] << "\t" << (integrator->GetSolver()->GetAcc())[1] << "\t" << (integrator->GetSolver()->GetAccVar())[0] << "\t" << (integrator->GetSolver()->GetAccVar())[1] << endl; - } - else{ - //if(currentFSIIter == config->GetStartTime()) historyFile << "\"FSI Iteration\"" << "\t" << "\"Displacement 1\"" << "\t" << "\"Displacement 1\"" << endl; - historyFile << currentFSIIter << "\t" << (integrator->GetSolver()->GetDisp())[0] << "\t" << (integrator->GetSolver()->GetDisp())[1] << endl; - } - } - } - - - if(config->GetUnsteady() == "YES"){ - if ( ExtIter%DeltaIter == 0 || ExtIter == NbExtIter ){ - restartFileName = config->GetRestartFile(); - restartFile.open(restartFileName.c_str(), ios::out); - if (structure->GetnDof() == 1){ - restartFile << "\"Time\"" << "\t" << "\"Displacement\"" << "\t" << "\"Velocity\"" << "\t" << "\"Acceleration\"" << "\t" << "\"Acceleration variable\"" << endl; - restartFile << currentTime-DeltaT << "\t" << (integrator->GetSolver()->GetDisp_n())[0] << "\t" << (integrator->GetSolver()->GetVel_n())[0] << "\t" << (integrator->GetSolver()->GetAcc_n())[0] << "\t" << (integrator->GetSolver()->GetAccVar_n())[0] << endl; - restartFile << currentTime << "\t" << (integrator->GetSolver()->GetDisp())[0] << "\t" << (integrator->GetSolver()->GetVel())[0] << "\t" << (integrator->GetSolver()->GetAcc())[0] << "\t" << (integrator->GetSolver()->GetAccVar())[0] << endl; - } - else if (structure->GetnDof() == 2){ - restartFile << "\"Time\"" << "\t" << "\"Displacement 1\"" << "\t" << "\"Displacement 2\"" << "\t" << "\"Velocity 1\"" << "\t" << "\"Velocity 2\"" << "\t" << "\"Acceleration 1\"" << "\t" << "\"Acceleration 2\"" << "\t" << "\"Acceleration variable 1\"" << "\t" << "\"Acceleration variable 2\"" << endl; - restartFile << currentTime-DeltaT << "\t" << (integrator->GetSolver()->GetDisp_n())[0] << "\t" << (integrator->GetSolver()->GetDisp_n())[1] << "\t" << (integrator->GetSolver()->GetVel_n())[0] << "\t" << (integrator->GetSolver()->GetVel_n())[1] << "\t" << (integrator->GetSolver()->GetAcc_n())[0] << "\t" << (integrator->GetSolver()->GetAcc_n())[1] << "\t" << (integrator->GetSolver()->GetAccVar_n())[0] << "\t" << (integrator->GetSolver()->GetAccVar_n())[1] << endl; - restartFile << currentTime << "\t" << (integrator->GetSolver()->GetDisp())[0] << "\t" << (integrator->GetSolver()->GetDisp())[1] << "\t" << (integrator->GetSolver()->GetVel())[0] << "\t" << (integrator->GetSolver()->GetVel())[1] << "\t" << (integrator->GetSolver()->GetAcc())[0] << "\t" << (integrator->GetSolver()->GetAcc())[1] << "\t" << (integrator->GetSolver()->GetAccVar())[0] << "\t" << (integrator->GetSolver()->GetAccVar())[1] << endl; - } - restartFile.close(); - } - } - -} - -void NativeSolidSolver::saveSolution(){ - - int rank = MASTER_NODE; - double DeltaT; - //string restartFileName; - //unsigned long DeltaIter = config->GetDeltaIterWrite(); - unsigned long ExtIter = integrator->GetExtIter(); - - DeltaT = config->GetDeltaT(); - - #ifdef HAVE_MPI - MPI_Comm_rank(MPI_COMM_WORLD, &rank); - #endif - - if(rank == MASTER_NODE){ - if(structure->GetnDof() == 1){ - if(config->GetUnsteady() == "YES"){ - historyFile << ExtIter << "\t" << (integrator->GetSolver()->GetDisp())[0] << "\t" << (integrator->GetSolver()->GetVel())[0] << "\t" << (integrator->GetSolver()->GetAcc())[0] << endl; - } - else{ - historyFile << ExtIter << "\t" << (integrator->GetSolver()->GetDisp())[0] << endl; - } - } - else if(structure->GetnDof() == 2){ - if(config->GetUnsteady() == "YES"){ - historyFile << ExtIter << "\t" << (integrator->GetSolver()->GetDisp())[0] << "\t" << (integrator->GetSolver()->GetDisp())[1] << "\t" << (integrator->GetSolver()->GetVel())[0] << "\t" << (integrator->GetSolver()->GetVel())[1] << "\t" << (integrator->GetSolver()->GetAcc())[0] << "\t" << (integrator->GetSolver()->GetAcc())[1] << endl; - } - else{ - historyFile << ExtIter << "\t" << (integrator->GetSolver()->GetDisp())[0] << "\t" << (integrator->GetSolver()->GetDisp())[1] << endl; - } - } - } -} - -void NativeSolidSolver::updateSolution(){ - - if(config->GetUnsteady() == "YES"){ - integrator->UpdateSolution(); - geometry->UpdateGeometry(); - structure->SetCenterOfRotation_n_X(structure->GetCenterOfRotation_x()); - structure->SetCenterOfRotation_n_Y(structure->GetCenterOfRotation_y()); - structure->SetCenterOfRotation_n_Z(structure->GetCenterOfRotation_z()); - } - else - q_uM1 = (integrator->GetSolver()->GetDisp()); -} - -/* - * void NativeSolidSolver::updateGeometry(){ - - if(config->GetUnsteady() == "YES"){ - geometry->UpdateGeometry(); - structure->SetCenterOfRotation_n_X(structure->GetCenterOfRotation_x()); - structure->SetCenterOfRotation_n_Y(structure->GetCenterOfRotation_y()); - structure->SetCenterOfRotation_n_Z(structure->GetCenterOfRotation_z()); - } -}*/ - -/* -void NativeSolidSolver::displacementPredictor(){ - - mapRigidBodyMotion(true, false); - -}*/ - -unsigned short NativeSolidSolver::getFSIMarkerID(){ - - unsigned short iMarker(0), IDtoSend; - - - while(iMarker < geometry->GetnMarkers()){ - if(geometry->GetMarkersMoving(iMarker)){ - IDtoSend = iMarker; - break; - } - iMarker++; - } - return IDtoSend; -} - -unsigned long NativeSolidSolver::getNumberOfSolidInterfaceNodes(unsigned short iMarker){ - return nSolidInterfaceVertex; -} - -unsigned int NativeSolidSolver::getInterfaceNodeGlobalIndex(unsigned short iMarker, unsigned short iVertex){ - - unsigned long iPoint; - - iPoint = geometry->vertex[iMarker][iVertex]; - - return iPoint; -} - -double NativeSolidSolver::getInterfaceNodePosX(unsigned short iMarker, unsigned short iVertex){ - - unsigned long iPoint; - double *Coord; - - iPoint = geometry->vertex[iMarker][iVertex]; - Coord = geometry->node[iPoint]->GetCoord(); - - return Coord[0]; -} - -double NativeSolidSolver::getInterfaceNodePosY(unsigned short iMarker, unsigned short iVertex){ - - unsigned long iPoint; - double *Coord; - - iPoint = geometry->vertex[iMarker][iVertex]; - Coord = geometry->node[iPoint]->GetCoord(); - - return Coord[1]; -} - -double NativeSolidSolver::getInterfaceNodePosZ(unsigned short iMarker, unsigned short iVertex){ - - unsigned long iPoint; - double *Coord; - - iPoint = geometry->vertex[iMarker][iVertex]; - Coord = geometry->node[iPoint]->GetCoord(); - - return 0.0; //3D is not really implemented in this solver... -} - -double NativeSolidSolver::getInterfaceNodePosX0(unsigned short iMarker, unsigned short iVertex){ - - unsigned long iPoint; - double *Coord; - - iPoint = geometry->vertex[iMarker][iVertex]; - Coord = geometry->node[iPoint]->GetCoord0(); - - return Coord[0]; -} - -double NativeSolidSolver::getInterfaceNodePosY0(unsigned short iMarker, unsigned short iVertex){ - - unsigned long iPoint; - double *Coord; - - iPoint = geometry->vertex[iMarker][iVertex]; - Coord = geometry->node[iPoint]->GetCoord0(); - - return Coord[1]; -} - -double NativeSolidSolver::getInterfaceNodePosZ0(unsigned short iMarker, unsigned short iVertex){ - - unsigned long iPoint; - double *Coord; - - iPoint = geometry->vertex[iMarker][iVertex]; - Coord = geometry->node[iPoint]->GetCoord0(); - - return 0.0; -} - -double NativeSolidSolver::getInterfaceNodeDispX(unsigned short iMarker, unsigned short iVertex){ - - unsigned long iPoint; - double *Coord, *Coord0; - - iPoint = geometry->vertex[iMarker][iVertex]; - Coord = geometry->node[iPoint]->GetCoord(); - Coord0 = geometry->node[iPoint]->GetCoord0(); - - return Coord[0] - Coord0[0]; -} - -double NativeSolidSolver::getInterfaceNodeDispY(unsigned short iMarker, unsigned short iVertex){ - - unsigned long iPoint; - double *Coord, *Coord0; - - iPoint = geometry->vertex[iMarker][iVertex]; - Coord = geometry->node[iPoint]->GetCoord(); - Coord0 = geometry->node[iPoint]->GetCoord0(); - return Coord[1] - Coord0[1]; -} - -double NativeSolidSolver::getInterfaceNodeDispZ(unsigned short iMarker, unsigned short iVertex){ - - unsigned long iPoint; - double *Coord, *Coord0; - - iPoint = geometry->vertex[iMarker][iVertex]; - Coord = geometry->node[iPoint]->GetCoord(); - Coord0 = geometry->node[iPoint]->GetCoord0(); - - return Coord[2] - Coord0[2]; -} - -double NativeSolidSolver::getInterfaceNodeVelX(unsigned short iMarker, unsigned short iVertex){ - unsigned long iPoint; - double *Vel; - - iPoint = geometry->vertex[iMarker][iVertex]; - Vel = geometry->node[iPoint]->GetVel(); - - return Vel[0]; -} - -double NativeSolidSolver::getInterfaceNodeVelY(unsigned short iMarker, unsigned short iVertex){ - unsigned long iPoint; - double *Vel; - - iPoint = geometry->vertex[iMarker][iVertex]; - Vel = geometry->node[iPoint]->GetVel(); - - return Vel[1]; -} - -double NativeSolidSolver::getInterfaceNodeVelZ(unsigned short iMarker, unsigned short iVertex){ - unsigned long iPoint; - double *Vel; - - iPoint = geometry->vertex[iMarker][iVertex]; - Vel = geometry->node[iPoint]->GetVel(); - - return Vel[2]; -} - -double NativeSolidSolver::getInterfaceNodeVelXNm1(unsigned short iMarker, unsigned short iVertex){ - unsigned long iPoint; - double *Vel_n; - - iPoint = geometry->vertex[iMarker][iVertex]; - Vel_n = geometry->node[iPoint]->GetVel_n(); - - return Vel_n[0]; -} - -double NativeSolidSolver::getInterfaceNodeVelYNm1(unsigned short iMarker, unsigned short iVertex){ - unsigned long iPoint; - double *Vel_n; - - iPoint = geometry->vertex[iMarker][iVertex]; - Vel_n = geometry->node[iPoint]->GetVel_n(); - - return Vel_n[1]; -} - -double NativeSolidSolver::getInterfaceNodeVelZNm1(unsigned short iMarker, unsigned short iVertex){ - unsigned long iPoint; - double *Vel_n; - - iPoint = geometry->vertex[iMarker][iVertex]; - Vel_n = geometry->node[iPoint]->GetVel_n(); - - return Vel_n[2]; -} - -double NativeSolidSolver::getRotationCenterPosX(){ - - return structure->GetCenterOfRotation_x(); - -} - -double NativeSolidSolver::getRotationCenterPosY(){ - - return structure->GetCenterOfRotation_y(); -} - -double NativeSolidSolver::getRotationCenterPosZ(){ - - return structure->GetCenterOfRotation_z(); -} - -void NativeSolidSolver::setGeneralisedForce(){ - - unsigned short iVertex, iMarker; - unsigned long iPoint; - double ForceX(0.0), ForceY(0.0), ForceZ(0.0); - double* force; - - iMarker = getFSIMarkerID(); - - for(iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++){ - iPoint = geometry->vertex[iMarker][iVertex]; - force = geometry->node[iPoint]->GetForce(); - ForceX += force[0]; - ForceY += force[1]; - ForceZ += force[2]; - } - - if(config->GetStructType() == "SPRING_HOR"){ - (integrator->GetSolver()->GetLoads())[0] = ForceX; - } - else if(config->GetStructType() == "SPRING_VER"){ - (integrator->GetSolver()->GetLoads())[0] = ForceY; - } - else if(config->GetStructType() == "AIRFOIL"){ - (integrator->GetSolver()->GetLoads())[0] = -ForceY; - } - else{ - cerr << "Wrong structural type for applying global fluild loads !" << endl; - throw(-1); - } - -} - -void NativeSolidSolver::setGeneralisedForce(double Fx, double Fy){ - - if(config->GetStructType() == "SPRING_HOR"){ - (integrator->GetSolver()->GetLoads())[0] = Fx; - } - else if(config->GetStructType() == "SPRING_VER"){ - (integrator->GetSolver()->GetLoads())[0] = Fy; - } - else if(config->GetStructType() == "AIRFOIL"){ - (integrator->GetSolver()->GetLoads())[0] = -Fy; - } - else{ - cerr << "Wrong structural type for applying global fluild loads !" << endl; - throw(-1); - } -} - -void NativeSolidSolver::setGeneralisedMoment(){ - - unsigned short iVertex, iMarker; - unsigned long iPoint; - double Moment(0.0), CenterX, CenterY, CenterZ; - double* Force; - double* Coord; - - iMarker = getFSIMarkerID(); - - for(iVertex = 0; iVertex < geometry->nVertex[iMarker]; iVertex++){ - iPoint = geometry->vertex[iMarker][iVertex]; - Force = geometry->node[iPoint]->GetForce(); - Coord = geometry->node[iPoint]->GetCoord(); - CenterX = getRotationCenterPosX(); - CenterY = getRotationCenterPosY(); - Moment += (Force[1]*(Coord[0]-CenterX) - Force[0]*(Coord[1]-CenterY)); - } - - if(config->GetStructType() == "AIRFOIL"){ - (integrator->GetSolver()->GetLoads())[1] = -Moment; - } - else if(config->GetStructType() == "SPRING_VER"){} - else if(config->GetStructType() == "SPRING_HOR"){} - else{ - cerr << "Wrong structural type for applying global fluild loads !" << endl; - throw(-1); - } - -} - -void NativeSolidSolver::setGeneralisedMoment(double M){ - - if(config->GetStructType() == "AIRFOIL"){ - (integrator->GetSolver()->GetLoads())[1] = -M; - } - else if(config->GetStructType() == "SPRING_VER"){} - else if(config->GetStructType() == "SPRING_HOR"){} - else{ - cerr << "Wrong structural type for applying global fluild loads !" << endl; - throw(-1); - } -} - -void NativeSolidSolver::applyload(unsigned short iVertex, double Fx, double Fy, double Fz){ - - unsigned short iMarker; - unsigned long iPoint; - double Force[3]; - - Force[0] = Fx; - Force[1] = Fy; - Force[2] = Fz; - - iMarker = getFSIMarkerID(); - iPoint = geometry->vertex[iMarker][iVertex]; - geometry->node[iPoint]->SetForce(Force); - -} diff --git a/api/NativeSolid_API.h b/api/NativeSolid_API.h deleted file mode 100644 index e67a995..0000000 --- a/api/NativeSolid_API.h +++ /dev/null @@ -1,75 +0,0 @@ -#pragma once - -#include -#include "../include/config.h" -#include "../include/structure.h" -#include "../include/integration.h" -#include "../include/output.h" -#include "../include/MatVec.h" -#include "../include/geometry.h" -#include -#include -#include - - -class NativeSolidSolver{ - -protected: - std::string confFile; - Config* config; - Geometry* geometry; - Structure* structure; - Integration* integrator; - Output* output; - std::ofstream historyFile; - std::ofstream restartFile; - CVector q_uM1; //The displacement at the previous FSI iteration - double omega; - unsigned long nSolidInterfaceVertex; - double varCoordNorm; - - -public: - /* NEW generation */ - NativeSolidSolver(std::string str, bool FSIComp); - ~NativeSolidSolver(); - void exit(); - //double getVarCoordNorm() const; - void preprocessIteration(unsigned long ExtIter); - void timeIteration(double t0, double tf); - //void mapRigidBodyMotion(bool predicition, bool initialize); - void computeInterfacePosVel(bool initialize); - void setInitialDisplacements(); - //void staticComputation(); - void writeSolution(double currentTime, double lastTime, double currentFSIIter, unsigned long ExtIter, unsigned long NbExtIter); - void saveSolution(); - void updateSolution(); - //void updateGeometry(); - //void displacementPredictor(); - unsigned short getFSIMarkerID(); - unsigned long getNumberOfSolidInterfaceNodes(unsigned short iMarker); - unsigned int getInterfaceNodeGlobalIndex(unsigned short iMarker, unsigned short iVertex); - double getInterfaceNodePosX(unsigned short iMarker, unsigned short iVertex); - double getInterfaceNodePosY(unsigned short iMarker, unsigned short iVertex); - double getInterfaceNodePosZ(unsigned short iMarker, unsigned short iVertex); - double getInterfaceNodePosX0(unsigned short iMarker, unsigned short iVertex); - double getInterfaceNodePosY0(unsigned short iMarker, unsigned short iVertex); - double getInterfaceNodePosZ0(unsigned short iMarker, unsigned short iVertex); - double getInterfaceNodeDispX(unsigned short iMarker, unsigned short iVertex); - double getInterfaceNodeDispY(unsigned short iMarker, unsigned short iVertex); - double getInterfaceNodeDispZ(unsigned short iMarker, unsigned short iVertex); - double getInterfaceNodeVelX(unsigned short iMarker, unsigned short iVertex); - double getInterfaceNodeVelY(unsigned short iMarker, unsigned short iVertex); - double getInterfaceNodeVelZ(unsigned short iMarker, unsigned short iVertex); - double getInterfaceNodeVelXNm1(unsigned short iMarker, unsigned short iVertex); - double getInterfaceNodeVelYNm1(unsigned short iMarker, unsigned short iVertex); - double getInterfaceNodeVelZNm1(unsigned short iMarker, unsigned short iVertex); - double getRotationCenterPosX(); - double getRotationCenterPosY(); - double getRotationCenterPosZ(); - void setGeneralisedForce(); - void setGeneralisedForce(double Fx, double Fy); - void setGeneralisedMoment(); - void setGeneralisedMoment(double M); - void applyload(unsigned short iVertex, double Fx, double Fy, double Fz); -}; diff --git a/include/MatVec.h b/include/MatVec.h deleted file mode 100644 index 364ffe9..0000000 --- a/include/MatVec.h +++ /dev/null @@ -1,106 +0,0 @@ -#pragma once - -#include "cblas.h" -#include "lapacke.h" - -#define ROW_MAJ 0 -#define COL_MAJ 1 - -class CVector { - -protected: - unsigned long nElm; - double* vec_val; - -public: - CVector(void); - CVector(const unsigned long & size, const double & val = 0.0); - CVector(const CVector & u); - ~CVector(); - unsigned long GetSize() const; - double* GetVec() const; - void print() const; - void Initialize(const unsigned long & size, const double & val = 0.0); - void SetAllValues(const double & val); - CVector & operator=(const CVector & u); - CVector & operator+=(const CVector & u); - CVector & operator-=(const CVector & u); - CVector & operator*=(const double & val); - CVector & operator/=(const double & val); - double & operator[](const unsigned long & i) const; - double norm() const; - double dotProd(const CVector & v) const; - void Reset(); -}; - -class CMatrix { - -protected: - unsigned long nEq; - unsigned long nVar; - CVector mat_val; - -public: - CMatrix(void); - CMatrix(const unsigned long & val_nEq, const unsigned long & val_nVar, const double & val = 0.0); - CMatrix(const CMatrix & A); - ~CMatrix(); - void Initialize(const unsigned long & val_nEq, const unsigned long & val_nVar, const double & val = 0.0); - void print() const; - CMatrix & operator=(const CMatrix & a); - CMatrix & operator+=(const CMatrix & a); - CMatrix & operator-=(const CMatrix & a); - CMatrix & operator*=(const double & val); - CMatrix & operator/=(const double & val); - unsigned long GetnEq() const; - unsigned long GetnVar() const; - double* GetMat() const; - CVector GetCVec() const; - void SetElm(const int & i, const int & j, double val); - double GetElm(const int & i, const int & j) const; - double DiagProduct() const; - double ComputeDet() const; - void Reset(); -}; - -CVector MatVecProd(const CMatrix & A, const CVector & b); -CVector ScalVecProd(const double & scal, const CVector & b); -CMatrix ScalMatProd(const double & scal, const CMatrix & A); -int SolveSys(const CMatrix & A, CVector & b); - -void MatrixToVec(int order, double** matrix, double* vecteur, int Nrow, int Ncol, int sizeVec); -void VecToMatrix(int order, double** matrix, double* vecteur, int Nrow, int Ncol, int sizeVec); - -/*OPERATOR +*/ -CVector operator+(CVector & vecA, CVector & vecB); -CVector operator+(const CVector & vecA, const CVector & vecB); -CVector operator+(CVector & vecA, const CVector & vecB); -CVector operator+(const CVector & vecA, CVector & vecB); -CMatrix operator+(const CMatrix & matA, const CMatrix & matB); - - -/*OPERATOR -*/ -CVector operator-(CVector & vecA, CVector & vecB); -CVector operator-(const CVector & vecA, CVector & vecB); -CVector operator-(CVector & vecA, const CVector & vecB); -CVector operator-(const CVector & vecA, const CVector & vecB); -CMatrix operator-(const CMatrix & matA, const CMatrix & matB); - -/*OPERATOR * */ -CVector operator*(CVector & vecA, double & val_mult); -CVector operator*(CVector & vecA, const double & val_mult); -CVector operator*(const CVector & vecA, double & val_mult); -CVector operator*(const CVector & vecA, const double & val_mult); -CVector operator*(double & val_mult, CVector & vecA); -CVector operator*(const double & val_mult, CVector & vecA); -CVector operator*(double & val_mult, const CVector & vecA); -CVector operator*(const double & val_mult, const CVector & vecA); -CMatrix operator*(CMatrix matA, double & val); -CMatrix operator*(double & val, CMatrix matA); - - -/*OPERATOR / */ -CVector operator/(CVector & vecA, double & val_mult); -CVector operator/(CVector & vecA, const double & val_mult); -CVector operator/(const CVector & vecA, double & val_mult); -CVector operator/(const CVector & vecA, const double & val_mult); diff --git a/include/config.h b/include/config.h deleted file mode 100644 index 087f6b8..0000000 --- a/include/config.h +++ /dev/null @@ -1,46 +0,0 @@ -#pragma once - -#include - - -class Config{ - -public: - Config(std::string filename); - virtual ~Config(); - virtual Config* GetAddress(); - virtual void ReadConfig(); - virtual std::string GetMeshFile(); - virtual std::string GetUnsteady(); - virtual std::string GetCSDSolver(); - virtual std::string GetStructType(); - virtual std::string GetLinearize(); - virtual std::string GetIntegrationAlgo(); - virtual std::string GetRestartSol(); - virtual std::string GetRestartFile(); - virtual std::string GetMovingMarker(); - virtual double GetSpringStiffness(); - virtual double GetSpringMass(); - virtual double GetInertiaCG(); - virtual double GetInertiaFlexural(); - virtual double GetSpringDamping(); - virtual double GetTorsionalStiffness(); - virtual double GetTorsionalDamping(); - virtual double GetCord(); - virtual double GetFlexuralAxis(); - virtual double GetGravityCenter(); - virtual double GetInitialDisp(); - virtual double GetInitialAngle(); - virtual double GetStartTime(); - virtual double GetDeltaT(); - virtual unsigned long GetDeltaIterWrite(); - virtual double GetStopTime(); - virtual double GetRho(); - -protected: - std::string ConfigFileName; - std::string MESH_FILE, UNSTEADY_SIMULATION, CSD_SOLVER, STRUCT_TYPE, LINEARIZE, INTEGRATION_ALGO, RESTART_SOL, RESTART_FILE, MOVING_MARKER; - double SPRING_STIFFNESS, SPRING_MASS, INERTIA_CG, INERTIA_FLEXURAL, SPRING_DAMPING, TORSIONAL_STIFFNESS, TORSIONAL_DAMPING, CORD, FLEXURAL_AXIS, GRAVITY_CENTER, INITIAL_DISP, INITIAL_ANGLE, START_TIME, DELTA_T, STOP_TIME, RHO; - unsigned long DELTAITERWRITE; - -}; diff --git a/include/geometry.h b/include/geometry.h deleted file mode 100644 index 2b255bd..0000000 --- a/include/geometry.h +++ /dev/null @@ -1,94 +0,0 @@ -#pragma once - -#include "config.h" -#include - -class Point{ - -protected: - double* Coord0; - double* Coord; - double* Coord_n; - double* Vel; - double* Vel_n; - double* Force; - -public: - Point(); - ~Point(); - double* GetCoord0() const; - double* GetCoord() const; - double* GetCoord_n() const; - double* GetVel() const; - double* GetVel_n() const; - double* GetForce() const; - void PrintCoord() const; - void SetCoord0(double* newCoord); - void SetCoord(double* newCoord); - void SetCoord_n(double* newCoord); - void SetVel(double* newVel); - void SetVel_n(double* newVel_n); - void SetForce(double* newForce); - void UpdateCoord(); - void UpdateVel(); -}; - -/*class Vertex{ - -protected: - double* Coord; - -public: - Vertex(); - ~Vertex(); - double* GetCoord() const; - void SetCoord(double* newCoord); - -}*/ - -/*class Element{ - -protected: - Point** pointElem; - unsigned short nPointElem; - -public: - Element(); - ~Element(); -}*/ - -/*class Marker{ - -protected: - std::string MarkerTag; - unsigned long nElem; - -public: - Marker(std::string Tag, unsigned long nElem_value); - ~Marker(); - Vertex** vertex; - Element** element; - -}*/ - -class Geometry{ - -protected: - unsigned short nDim; - unsigned long nElem, nPoint, nMarkers; - unsigned long* nElemMarker; - -public: - Geometry(Config* config); - ~Geometry(); - unsigned long GetnMarkers() const; - bool GetMarkersMoving(unsigned long iMarker) const; - void UpdateGeometry(); - unsigned long** vertex; //vertex[iMarker][iPoint] - unsigned long* nVertex; //nVertex[iMarker] - bool* markersMoving; - Point** node; //node[iPoint] returns a pointeur to a Point object so e.g. node[iPoint]->PrintCoord(); -}; - -bool isInVec(std::vector const& inputVector, int dummyInt); -void vecCopy(std::vector const& sourceVector, unsigned long* destinationTab); diff --git a/include/integration.h b/include/integration.h deleted file mode 100644 index 0d243d1..0000000 --- a/include/integration.h +++ /dev/null @@ -1,33 +0,0 @@ -#pragma once - -//#include "MatVec.h" -#include "config.h" -#include "structure.h" -#include "solver.h" - -#define PI 3.14159265 - - -class Integration{ - -protected: - double totTime; - double deltaT; - unsigned long ExtIter; - std::string algo; - - Solver* solver; - -public: - Integration(Config* config, Structure* structure); - ~Integration(); - Solver* GetSolver(); - double GettotTime(); - double GetdeltaT(); - void SetExtIter(unsigned long val_ExtIter); - unsigned long GetExtIter(); - void SetInitialConditions(Config* config, Structure* structure); - void TemporalIteration(double& t0, double& tf, Structure *structure); - //void StaticIteration(Config *config, Structure *structure); - void UpdateSolution(); -}; diff --git a/include/output.h b/include/output.h deleted file mode 100644 index 44ba131..0000000 --- a/include/output.h +++ /dev/null @@ -1,19 +0,0 @@ -#pragma once - -#include - -#include "MatVec.h" -#include "config.h" -#include "structure.h" -#include "integration.h" - -class Output{ - -public: - Output(void); - ~Output(); - //void WriteHistory(Integration* solver, Structure* structure, std::ofstream* outputfile, const double & time); - //void WriteRestart(Integration* solver, Structure* structure); - //void WriteStaticSolution(Config* config, Integration* solver, Structure* structure, std::ofstream* outputfile); - //void WriteRestart(); -}; diff --git a/include/solver.h b/include/solver.h deleted file mode 100644 index f806e52..0000000 --- a/include/solver.h +++ /dev/null @@ -1,85 +0,0 @@ -#pragma once - -#include "MatVec.h" -#include "structure.h" -#include "config.h" - -class Solver{ - -protected: - CVector q; - CVector qdot; - CVector qddot; - CVector q_n; - CVector qdot_n; - CVector qddot_n; - CVector Loads; - CVector Loads_n; - CVector a; - CVector a_n; - bool linear; - -public: - Solver(unsigned int nDof, bool bool_linear); - virtual ~Solver(); - virtual void Iterate(double& t0, double& tf, Structure* structure); - virtual CVector & GetDisp(); - virtual CVector & GetVel(); - virtual CVector & GetAcc(); - virtual CVector & GetDisp_n(); - virtual CVector & GetVel_n(); - virtual CVector & GetAcc_n(); - virtual CVector & GetLoads(); - virtual CVector & GetAccVar(); - virtual CVector & GetAccVar_n(); - virtual void ResetSolution(); - virtual void SaveToThePast(); - virtual void SetInitialState(Config *config, Structure* structure); - -}; - -class AlphaGenSolver : public Solver { - -protected: - double beta; - double gamma; - double alpha_m; - double alpha_f; - double rho; - double gammaPrime; - double betaPrime; - -public: - AlphaGenSolver(unsigned int nDof, double val_rho, bool bool_linear); - ~AlphaGenSolver(); - CVector & GetAccVar(); - CVector & GetAccVar_n(); - virtual void Iterate(double& t0, double& tf, Structure* structure); - void ComputeRHS(Structure* structure, CVector &RHS); - void ComputeResidual(Structure* structure, CVector &res); - void ComputeTangentOperator(Structure* structure, CMatrix & St); - void ResetSolution(); - void SaveToThePast(); - virtual void SetInitialState(Config *config, Structure* structure); - -}; - - -class RK4Solver : public Solver { - -protected: - unsigned int size; - double lastTime; - double currentTime; - -public: - RK4Solver(unsigned nDof, bool bool_linear); - ~RK4Solver(); - virtual void Iterate(double &t0, double &tf, Structure* structure); - void EvaluateStateDerivative(double tCurrent, CVector& state, CVector& stateDerivative, Structure* structure); - void interpLoads(double& tCurrent, CVector& val_loads); - virtual void SetInitialState(Config* config, Structure* structure); - CVector SetState(); - CVector SetState_n(); - -}; diff --git a/include/structure.h b/include/structure.h deleted file mode 100644 index 5f473f5..0000000 --- a/include/structure.h +++ /dev/null @@ -1,52 +0,0 @@ -#pragma once - -#include "MatVec.h" -#include "config.h" - - -class Structure{ - -protected: - unsigned int nDof; - double m; - double Kh; - double Ka; - double Ch; - double Ca; - double c; - double b; - double xf; - double xCG; - double S; - double ICG; - double If; - double centerOfRotation[3]; - double centerOfRotation_n[3]; - - -public: - Structure(Config* config); - virtual ~Structure(); - void SetCenterOfRotation_X(double coord_x); - void SetCenterOfRotation_Y(double coord_y); - void SetCenterOfRotation_Z(double coord_z); - double GetCenterOfRotation_x() const; - double GetCenterOfRotation_y() const; - double GetCenterOfRotation_z() const; - const double* GetCenterOfRotation() const; - void SetCenterOfRotation_n_X(double coord_x); - void SetCenterOfRotation_n_Y(double coord_y); - void SetCenterOfRotation_n_Z(double coord_z); - double GetCenterOfRotation_n_x() const; - double GetCenterOfRotation_n_y() const; - double GetCenterOfRotation_n_z() const; - unsigned int GetnDof() const; - double Get_m() const; - double Get_Kh() const; - double Get_Ka() const; - double Get_Ch() const; - double Get_Ca() const; - double Get_S() const; - double Get_If() const; - -}; diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt deleted file mode 100644 index 8c16a27..0000000 --- a/src/CMakeLists.txt +++ /dev/null @@ -1,9 +0,0 @@ -# Compilation of the source code -# Generation of a shared librairy libNative.so - -ADD_LIBRARY(Native SHARED ../include/config.h ../include/geometry.h ../include/integration.h ../include/MatVec.h ../include/output.h ../include/structure.h ../include/solver.h config.cpp geometry.cpp integration.cpp MatVec.cpp output.cpp structure.cpp solver.cpp ../api/NativeSolid_API.h ../api/NativeSolid_API.cpp) - -TARGET_LINK_LIBRARIES(Native ${LAPACKE_LIBRARIES} ${CBLAS_LIBRARIES} ${BLAS_LIBRARIES}) - -ADD_EXECUTABLE(TestCVector ../include/MatVec.h MatVec.cpp testcvector.cpp) -TARGET_LINK_LIBRARIES(TestCVector ${LAPACKE_LIBRARIES} ${CBLAS_LIBRARIES} ${BLAS_LIBRARIES}) diff --git a/src/MatVec.cpp b/src/MatVec.cpp deleted file mode 100644 index 2f614af..0000000 --- a/src/MatVec.cpp +++ /dev/null @@ -1,579 +0,0 @@ -#include "../include/MatVec.h" -#include -#include - -using namespace std; - -/*--- Definition of the CVector class ---*/ - -CVector::CVector(void):nElm(0){ - - vec_val = NULL; - -} - -CVector::CVector(const unsigned long & size, const double & val){ - - vec_val = NULL; - if(size <= 0){ - nElm = 0; - cerr << "CVector:CVector(const unsigned long &, const double): " << "invalid input : size = " << size << endl; - throw(-1); - } - else{ - nElm = size; - vec_val = new double[nElm]; - for(unsigned int i=0; i < nElm; i++) - vec_val[i] = val; - } -} - -CVector::CVector(const CVector & u){ - - nElm = u.nElm; - vec_val = NULL; - vec_val = new double[nElm]; - for(unsigned int i = 0; i < nElm; i++) - vec_val[i] = u.vec_val[i]; -} - -CVector::~CVector(){ - - if (vec_val != NULL) delete [] vec_val; -} - -void CVector::Initialize(const unsigned long &size, const double & val){ - - if(vec_val == NULL){ - if(size <= 0){ - nElm = 0; - cerr << "CVector::Initialize(const unsigned long &, const double &): " << "invalid number of element" << nElm << endl; - throw(-1); - } - else{ - nElm = size; - vec_val = new double[nElm]; - for(unsigned int ii=0; ii -#include -#include -#include -#include - -using namespace std; - -Config::Config(string filename):ConfigFileName(filename) -{ - -} - -Config::~Config() -{ - -} - -Config* Config::GetAddress() -{ - return this; -} - -void Config::ReadConfig() -{ - string delimiter = "="; - size_t pos; - string text_line, option; - ifstream InputFile; - InputFile.open(ConfigFileName.c_str(), ios::in); - if(InputFile.fail()){ - cerr << "Invalid configuration file name : " << ConfigFileName << endl; - throw(-1); - } - char caract; - while (getline(InputFile,text_line)) - { - caract = text_line[0]; - if (caract == '%') - { - } - else - { - pos = text_line.find(delimiter); - option = text_line.substr(0,pos); - text_line.erase(0,pos+delimiter.length()); - option.erase(remove(option.begin(), option.end(), ' '), option.end()); - text_line.erase(remove(text_line.begin(), text_line.end(),' '),text_line.end()); - if (option == "CSD_SOLVER") CSD_SOLVER = text_line; - else if (option == "MESH_FILE") MESH_FILE = text_line; - else if (option == "UNSTEADY_SIMULATION") UNSTEADY_SIMULATION = text_line; - else if (option == "STRUCT_TYPE") STRUCT_TYPE = text_line; - else if (option == "LINEARIZE") LINEARIZE = text_line; - else if (option == "INTEGRATION_ALGO") INTEGRATION_ALGO = text_line; - else if (option == "RESTART_SOL") RESTART_SOL = text_line; - else if (option == "RESTART_FILE") RESTART_FILE = text_line; - else if (option == "MOVING_MARKER") MOVING_MARKER = text_line; - else if (option == "SPRING_MASS") SPRING_MASS = atof(text_line.c_str()); - else if (option == "INERTIA_CG") INERTIA_CG = atof(text_line.c_str()); - else if (option == "INERTIA_FLEXURAL") INERTIA_FLEXURAL = atof(text_line.c_str()); - else if (option == "SPRING_STIFFNESS") SPRING_STIFFNESS = atof(text_line.c_str()); - else if (option == "SPRING_DAMPING") SPRING_DAMPING = atof(text_line.c_str()); - else if (option == "TORSIONAL_STIFFNESS") TORSIONAL_STIFFNESS = atof(text_line.c_str()); - else if (option == "TORSIONAL_DAMPING") TORSIONAL_DAMPING = atof(text_line.c_str()); - else if (option == "CORD") CORD = atof(text_line.c_str()); - else if (option == "FLEXURAL_AXIS") FLEXURAL_AXIS = atof(text_line.c_str()); - else if (option == "GRAVITY_CENTER") GRAVITY_CENTER = atof(text_line.c_str()); - else if (option == "INITIAL_DISP") INITIAL_DISP = atof(text_line.c_str()); - else if (option == "INITIAL_ANGLE") INITIAL_ANGLE = atof(text_line.c_str()); - else if (option == "START_TIME") START_TIME = atof(text_line.c_str()); - else if (option == "DELTA_T") DELTA_T = atof(text_line.c_str()); - else if (option == "DELTA_ITER_WRITE") DELTAITERWRITE = atol(text_line.c_str()); - else if (option == "STOP_TIME") STOP_TIME = atof(text_line.c_str()); - else if (option == "RHO") RHO = atof(text_line.c_str()); - else cout << "The option " + option + " is not recognized !" << endl; - } - } - InputFile.close(); - - if (CSD_SOLVER == "NATIVE") cout << "The Native solver has been chosen" << endl; - else cout << "Cannot run the solver with other value than NATIVE for CSD_SOLVER option !" << endl; - - if (UNSTEADY_SIMULATION == "YES") cout << "Dynamic structure computation" << endl; - else cout << "Static structure computation" << endl; - - if (STRUCT_TYPE == "SPRING_HOR" || STRUCT_TYPE == "SPRING_VER") cout << "Structural model is a plunging spring" << endl; - else if (STRUCT_TYPE == "AIRFOIL") cout << "Structural model is a pitching plunging airfoil" << endl; - else cout << "The specified structural model is not recognized or implemented yet !" << endl; -} - -std::string Config::GetMeshFile() -{ - return MESH_FILE; -} - -std::string Config::GetUnsteady() -{ - return UNSTEADY_SIMULATION; -} - -std::string Config::GetCSDSolver() -{ - return CSD_SOLVER; -} - -std::string Config::GetStructType() -{ - return STRUCT_TYPE; -} - -std::string Config::GetLinearize(){ - return LINEARIZE; -} - -std::string Config::GetIntegrationAlgo() -{ - return INTEGRATION_ALGO; -} - -std::string Config::GetRestartSol() -{ - return RESTART_SOL; -} - -std::string Config::GetRestartFile() -{ - return RESTART_FILE; -} - -std::string Config::GetMovingMarker() -{ - return MOVING_MARKER; -} - -double Config::GetStartTime(){ - return START_TIME; -} - -double Config::GetDeltaT(){ - return DELTA_T; -} - -unsigned long Config::GetDeltaIterWrite(){ - return DELTAITERWRITE; -} - -double Config::GetStopTime(){ - return STOP_TIME; -} - -double Config::GetSpringStiffness(){ - return SPRING_STIFFNESS; -} - -double Config::GetSpringMass(){ - return SPRING_MASS; -} - -double Config::GetInertiaCG(){ - return INERTIA_CG; -} - -double Config::GetInertiaFlexural(){ - return INERTIA_FLEXURAL; -} - -double Config::GetSpringDamping(){ - return SPRING_DAMPING; -} - -double Config::GetTorsionalStiffness(){ - return TORSIONAL_STIFFNESS; -} - -double Config::GetTorsionalDamping(){ - return TORSIONAL_DAMPING; -} - -double Config::GetCord(){ - return CORD; -} - -double Config::GetFlexuralAxis(){ - return FLEXURAL_AXIS; -} - -double Config::GetGravityCenter(){ - return GRAVITY_CENTER; -} - -double Config::GetInitialDisp(){ - return INITIAL_DISP; -} - -double Config::GetInitialAngle(){ - return INITIAL_ANGLE; -} - -double Config::GetRho(){ - return RHO; -} diff --git a/src/geometry.cpp b/src/geometry.cpp deleted file mode 100644 index 105d525..0000000 --- a/src/geometry.cpp +++ /dev/null @@ -1,331 +0,0 @@ -#include "../include/geometry.h" -#include "../include/config.h" -#include -#include -#include -#include -#include -#include - -using namespace std; - -Point::Point(){ - - Coord0 = new double[3]; - Coord0[0] = 0.0; - Coord0[1] = 0.0; - Coord0[2] = 0.0; - - Coord = new double[3]; - Coord[0] = 0.0; - Coord[1]= 0.0; - Coord[2] = 0.0; - - Coord_n = new double[3]; - Coord_n[0] = 0.0; - Coord_n[1] = 0.0; - Coord_n[2] = 0.0; - - Vel = new double[3]; - Vel[0] = 0.0; - Vel[1] = 0.0; - Vel[2] = 0.0; - - Vel_n = new double[3]; - Vel_n[0] = 0.0; - Vel_n[1] = 0.0; - Vel_n[2] = 0.0; - - Force = new double[3]; - Force[0] = 0.0; - Force[1] = 0.0; - Force[2] = 0.0; - -} - -Point::~Point(){ - delete [] Coord0; - Coord0 = NULL; - - delete [] Coord; - Coord = NULL; - - delete [] Coord_n; - Coord_n = NULL; - - delete [] Vel; - Vel = NULL; - - delete [] Vel_n; - Vel_n = NULL; - - delete [] Force; - Force = NULL; -} - -double* Point::GetCoord0() const{ - return Coord0; -} - -double* Point::GetCoord() const{ - return Coord; -} - -double* Point::GetCoord_n() const{ - return Coord_n; -} - -double* Point::GetVel() const{ - return Vel; -} - -double* Point::GetVel_n() const{ - return Vel_n; -} - -double* Point::GetForce() const{ - return Force; -} - -void Point::PrintCoord() const{ - cout << Coord[0] << " ; " << Coord[1] << " ; " << Coord[2] << endl; -} - -void Point::SetCoord0(double* newCoord){ - Coord0[0] = newCoord[0]; - Coord0[1] = newCoord[1]; - Coord0[2] = newCoord[2]; -} - -void Point::SetCoord(double* newCoord){ - Coord[0] = newCoord[0]; - Coord[1] = newCoord[1]; - Coord[2] = newCoord[2]; -} - -void Point::SetCoord_n(double* newCoord){ - Coord_n[0] = newCoord[0]; - Coord_n[1] = newCoord[1]; - Coord_n[2] = newCoord[2]; -} - -void Point::SetVel(double *newVel){ - Vel[0] = newVel[0]; - Vel[1] = newVel[1]; - Vel[2] = newVel[2]; -} - -void Point::SetVel_n(double *newVel_n){ - Vel_n[0] = newVel_n[0]; - Vel_n[1] = newVel_n[1]; - Vel_n[2] = newVel_n[2]; -} - -void Point::SetForce(double* newForce){ - Force[0] = newForce[0]; - Force[1] = newForce[1]; - Force[2] = newForce[2]; -} - -void Point::UpdateCoord(){ - Coord_n[0] = Coord[0]; - Coord_n[1] = Coord[1]; - Coord_n[2] = Coord[2]; -} - -void Point::UpdateVel(){ - Vel_n[0] = Vel[0]; - Vel_n[1] = Vel[1]; - Vel_n[2] = Vel[2]; -} - -Geometry::Geometry(Config* config){ - - string meshFileName, textLine, tampon; - ifstream meshFile; - string::size_type position; - double Coord[3]; - double* TempCoord; - unsigned long iMarker(0); - int elemType(0), dummyInt(0); - int iPoint; - - Coord[0] = 0.0; - Coord[1] = 0.0; - Coord[2] = 0.0; - - nDim = 0; - nElem = 0; - nPoint = 0; - nMarkers = 0; - //strcpy(cstr, config->GetMeshFile()); - - meshFileName = config->GetMeshFile(); - - /*--- Open the mesh file and check ---*/ - meshFile.open(meshFileName.c_str(), ios::in); - if (meshFile.fail()){ - cout << "Enable to open the mesh file " << meshFileName << endl; - exit(1); - } - - cout << "Mesh file " << meshFileName << " is open." << endl; - cout << "Constructing the mesh..." << endl; - while(getline(meshFile, textLine)){ - - position = textLine.find("NDIME=",0); - if (position != string::npos){ - textLine.erase(0,6); - nDim = atoi(textLine.c_str()); - cout << "Number of dimensions : " << nDim << endl; - } - - position = textLine.find("NELEM=",0); - if (position != string::npos){ - textLine.erase(0,6); - nElem = atoi(textLine.c_str()); - cout << "Number of elements : " << nElem << endl; - for(int iElem=0; iElem < nElem; iElem++){ - getline(meshFile, textLine); - } - } - - position = textLine.find("NPOIN=", 0); - if (position != string::npos){ - textLine.erase(0,6); - nPoint = atoi(textLine.c_str()); - cout << "Number of points : " << nPoint << endl; - node = new Point*[nPoint]; - //cout << "JE VAIS REMPLIR" << endl; - for(iPoint=0; iPoint < nPoint; iPoint++){ - getline(meshFile, textLine); - //if(iPoint == 1) cout << textLine << endl; - node[iPoint] = new Point(); - istringstream point_line(textLine); - point_line >> Coord[0]; - //if(iPoint == 1) cout << Coord[0] << endl; - point_line >> Coord[1]; - //if(iPoint == 1) cout << Coord[1] << endl; - if(nDim == 3) point_line >> Coord[2]; - node[iPoint]->SetCoord0(Coord); - node[iPoint]->SetCoord(Coord); - node[iPoint]->SetCoord_n(Coord); - TempCoord = node[iPoint]->GetCoord(); - //cout << iPoint << endl; - } - } - - position = textLine.find("NMARK=", 0); - if (position != string::npos){ - textLine.erase(0,6); - nMarkers = atoi(textLine.c_str()); - cout << "Number of markers : " << nMarkers << endl; - vertex = new unsigned long*[nMarkers]; - nVertex = new unsigned long[nMarkers]; - markersMoving = new bool[nMarkers]; - nElemMarker = new unsigned long[nMarkers]; - } - - position = textLine.find("MARKER_TAG=",0); - if (position != string::npos){ - textLine.erase(0,12); - cout << "Reading elements for marker : " << textLine << endl; - if (textLine == config->GetMovingMarker()){ - cout << "Marker " << textLine << " is a moving marker." << endl; - markersMoving[iMarker] = true; - } - else markersMoving[iMarker] = false; - } - - position = textLine.find("MARKER_ELEMS=",0); - if (position != string::npos){ - textLine.erase(0,13); - nElemMarker[iMarker] = atoi(textLine.c_str()); - cout << "Number of elements on the marker : " << nElemMarker[iMarker] << endl; - vector tempVertexMarker; - for(int iElem=0; iElem < nElemMarker[iMarker]; iElem++){ - getline(meshFile,textLine); - istringstream point_line(textLine); - point_line >> elemType; - if (elemType == 3){ - point_line >> dummyInt; - if (!isInVec(tempVertexMarker, dummyInt)){ - tempVertexMarker.push_back(dummyInt); - } - point_line >> dummyInt; - if (!isInVec(tempVertexMarker, dummyInt)){ - tempVertexMarker.push_back(dummyInt); - } - } - else{ - cout << "Elem type " << elemType << " is not recognized !" << endl; - exit(1); - } - } - nVertex[iMarker] = tempVertexMarker.size(); - vertex[iMarker] = new unsigned long[nVertex[iMarker]]; - vecCopy(tempVertexMarker,vertex[iMarker]); - iMarker++; - } - - } - - meshFile.close(); - - /*for (int iVertex = 0; iVertex < nVertex[0]; iVertex++){ - //cout << vertex[0][iVertex] << endl; - iPoint = vertex[0][iVertex]; - node[iPoint]->PrintCoord(); - }*/ - - - -} - -Geometry::~Geometry(){ - - for(int iPoint=0; iPoint < nPoint; iPoint++){ - delete node[iPoint]; - } - - for(int iMarker=0; iMarker < nMarkers; iMarker++){ - delete [] vertex[iMarker]; - } - - delete [] vertex; - - delete [] node; - - delete [] nElemMarker; - -} - -unsigned long Geometry::GetnMarkers() const{ - return nMarkers; -} - -bool Geometry::GetMarkersMoving(unsigned long iMarker) const{ - return markersMoving[iMarker]; -} - -void Geometry::UpdateGeometry(){ - - unsigned long iPoint; - - for(iPoint=0; iPoint < nPoint; iPoint++){ - node[iPoint]->UpdateCoord(); - node[iPoint]->UpdateVel(); - } -} - -bool isInVec(vector const& inputVector, int dummyInt){ - int i = 0; - while(i const& sourceVector, unsigned long* destinationTab){ - for(int i=0; i -#include -#include -#include -#include - -using namespace std; - -Integration::Integration(Config *config, Structure *structure){ - - solver = NULL; - bool linear; - - linear = (config->GetLinearize()) == "YES"; - - if(config->GetIntegrationAlgo() == "ALPHAGEN"){ - solver = new AlphaGenSolver(structure->GetnDof(), config->GetRho(), linear); - } - else if(config->GetIntegrationAlgo() == "RK4"){ - solver = new RK4Solver(structure->GetnDof(), linear); - } - else{ - - } -} - -Integration::~Integration(){ - - if (solver != NULL) delete solver; -} - -Solver* Integration::GetSolver(){ - - return solver; -} - -double Integration::GettotTime(){ - return totTime; -} - -double Integration::GetdeltaT(){ - return deltaT; -} - -void Integration::SetExtIter(unsigned long val_ExtIter){ - - ExtIter = val_ExtIter; -} - -unsigned long Integration::GetExtIter(){ - - return ExtIter; -} - -void Integration::SetInitialConditions(Config* config, Structure* structure){ - - ExtIter = 0; - solver->SetInitialState(config, structure); - structure->SetCenterOfRotation_Y((solver->GetDisp())[0]); - -} - -void Integration::TemporalIteration(double& t0, double& tf, Structure *structure){ - - int rank = 0; - -#ifdef HAVE_MPI - MPI_Comm_rank(MPI_COMM_WORLD, &rank); -#endif - - double currentTime(tf); - deltaT = tf - t0; - solver->Iterate(t0, tf, structure); - - if(rank == 0){ - if(structure->GetnDof() == 1){ - cout << " Time" << "\t" << "Displacement" << "\t" << "Velocity" << "\t" << "Acceleration" << "\t" << endl; - cout << " " << currentTime << "\t" << (solver->GetDisp())[0] << "\t" << (solver->GetVel())[0] << "\t" << (solver->GetAcc())[0] << endl; - } - else if(structure->GetnDof() == 2){ - cout << " Time" << "\t" << "Displacement 1" << "\t" << "Displacement 2" << "\t" << "Velocity 1" << "\t" << "Velocity 2" << "\t" << "Acceleration 1" << "\t" << "Acceleration 2" << endl; - cout << currentTime << "\t" << (solver->GetDisp())[0] << "\t" << (solver->GetDisp())[1] << "\t" << (solver->GetVel())[0] << "\t" << (solver->GetVel())[1] << "\t" << (solver->GetAcc())[0] << "\t" << (solver->GetAcc())[1] << endl; - } - } -} - -/* -void Integration::StaticIteration(Config *config, Structure *structure){ - q->Reset(); - *q += *Loads; - SolveSys(structure->GetK(),q); //q will be updated with the new solution... -} -*/ - -void Integration::UpdateSolution(){ - - solver->SaveToThePast(); - -} diff --git a/src/output.cpp b/src/output.cpp deleted file mode 100644 index da7bc8a..0000000 --- a/src/output.cpp +++ /dev/null @@ -1,63 +0,0 @@ -#include "../include/output.h" -#include - -using namespace std; - -Output::Output(){} - -Output::~Output(){} - -/* -void Output::WriteHistory(Integration* solver, Structure* structure, ofstream* outputfile, const double & time){ - - if(structure->GetnDof() == 1){ - if(time == 0){ - cout << "\"Time\"" << "\t" << "\"Displacement\"" << "\t" << "\"Velocity\"" << "\t" << "\"Acceleration\"" << "\t" << endl; - outputfile[0] << "\"Time\"" << "\t" << "\"Displacement\"" << "\t" << "\"Velocity\"" << "\t" << "\"Acceleration\"" << "\t" << endl; - } - cout << time << "\t" << (*(solver->GetDisp()))[0] << "\t" << (*(solver->GetVel()))[0] << "\t" << (*(solver->GetAcc()))[0] << endl; - outputfile[0] << time << "\t" << (*(solver->GetDisp()))[0] << "\t" << (*(solver->GetVel()))[0] << "\t" << (*(solver->GetAcc()))[0] << endl; - } - else if(structure->GetnDof() == 2){ - if(time == 0){ - cout << "\"Time\"" << "\t" << "\"Displacement 1\"" << "\t" << "\"Displacement 2\"" << "\t" << "\"Velocity 1\"" << "\t" << "\"Velocity 2\"" << "\t" << "\"Acceleration 1\"" << "\t" << "\"Acceleration 2\"" << endl; - outputfile[0] << "\"Time\"" << "\t" << "\"Displacement 1\"" << "\t" << "\"Displacement 2\"" << "\t" << "\"Velocity 1\"" << "\t" << "\"Velocity 2\"" << "\t" << "\"Acceleration 1\"" << "\t" << "\"Acceleration 2\"" << endl; - } - cout << time << "\t" << (*(solver->GetDisp()))[0] << "\t" << (*(solver->GetDisp()))[1] << "\t" << (*(solver->GetVel()))[0] << "\t" << (*(solver->GetVel()))[1] << "\t" << (*(solver->GetAcc()))[0] << "\t" << (*(solver->GetAcc()))[1] << endl; - outputfile[0] << time << "\t" << (*(solver->GetDisp()))[0] << "\t" << (*(solver->GetDisp()))[1] << "\t" << (*(solver->GetVel()))[0] << "\t" << (*(solver->GetVel()))[1] << "\t" << (*(solver->GetAcc()))[0] << "\t" << (*(solver->GetAcc()))[1] << endl; - } -} - -void Output::WriteRestart(Integration* solver, Structure* structure){ - ofstream RestartFile; - RestartFile.open("Nat_solution_restart.out", ios::out); - - if(structure->GetnDof() == 1){ - RestartFile << "\"Displacement\"" << "\t" << "\"Velocity\"" << "\t" << "\"Acceleration\"" << "\t" << "\"Acceleration variable\"" << endl; - RestartFile << (*(solver->GetDisp()))[0] << "\t" << (*(solver->GetVel()))[0] << "\t" << (*(solver->GetAcc()))[0] << "\t" << (*(solver->GetAccVar()))[0] << endl; - } - else if(structure->GetnDof() == 2){ - RestartFile << "\"Displacement 1\"" << "\t" << "\"Displacement 2\"" << "\t" << "\"Velocity 1\"" << "\t" << "\"Velocity 2\"" << "\t" << "\"Acceleration 1\"" << "\t" << "\"Acceleration 2\"" << "\t" << "\"Acceleration variable 1\"" << "\t" << "\"Acceleration variable 2\"" << endl; - RestartFile << time << "\t" << (*(solver->GetDisp()))[0] << "\t" << (*(solver->GetDisp()))[1] << "\t" << (*(solver->GetVel()))[0] << "\t" << (*(solver->GetVel()))[1] << "\t" << (*(solver->GetAcc()))[0] << "\t" << (*(solver->GetAcc()))[1]<< "\t" << (*(solver->GetAccVar()))[0] << "\t" << (*(solver->GetAccVar()))[1] << endl; - } -} - -void Output::WriteStaticSolution(Config* config, Integration* solver, Structure* structure, ofstream* outputfile){ - if(structure->GetnDof() == 1){ - cout << "Static displacement is : " << (*(solver->GetDisp()))[0] << " [m]" << endl; - cout << "Writing displacement into a solution file" << endl; - if(config->GetStructType() == "SPRING_HOR"){ - outputfile[0] << -1.000 << "\t" << (*(solver->GetDisp()))[0] << "\t" << 0.000; - } - else if(config->GetStructType() == "SPRING_VER"){ - outputfile[0] << -1.000 << "\t" << 0.000 << "\t" << (*(solver->GetDisp()))[0]; - } - } - else if(structure->GetnDof() == 2){ - cout << "Plunging displacement is :" << (*(solver->GetDisp()))[0] << " [m]" << endl; - cout << "Pitching rotation is :" << (*(solver->GetDisp()))[1] << " [rad]" << endl; - cout << "Writing displacement into a solution file" << endl; - outputfile[0] << -1.000 << "\t" << (*(solver->GetDisp()))[0] << "\t" << (*(solver->GetDisp()))[1]; - } -} -*/ diff --git a/src/solver.cpp b/src/solver.cpp deleted file mode 100644 index 6ae3af9..0000000 --- a/src/solver.cpp +++ /dev/null @@ -1,648 +0,0 @@ -#include "../include/solver.h" -#include "../include/structure.h" -#include "../include/MatVec.h" - -#include -#include -#include - -using namespace std; - -/* CLASS SOLVER*/ -Solver::Solver(unsigned int nDof, bool bool_linear){ - - q.Initialize(nDof, 0.0); - qdot.Initialize(nDof, 0.0); - qddot.Initialize(nDof, 0.0); - q_n.Initialize(nDof, 0.0); - qdot_n.Initialize(nDof, 0.0); - qddot_n.Initialize(nDof, 0.0); - Loads.Initialize(nDof, 0.0); - Loads_n.Initialize(nDof, 0.0); - - ResetSolution(); - Loads.Reset(); - Loads_n.Reset(); - - linear = bool_linear; - -} - -Solver::~Solver(){} - -void Solver::Iterate(double &t0, double &tf, Structure* structure){} - -CVector & Solver::GetDisp(){ - return q; -} - -CVector & Solver::GetVel(){ - return qdot; -} - -CVector & Solver::GetAcc(){ - return qddot; -} - -CVector & Solver::GetDisp_n(){ - return q_n; -} - -CVector & Solver::GetVel_n(){ - return qdot_n; -} - -CVector & Solver::GetAcc_n(){ - return qddot_n; -} - -CVector & Solver::GetLoads(){ - return Loads; -} - -CVector & Solver::GetAccVar(){ - return a; -} - -CVector & Solver::GetAccVar_n(){ - return a_n; -} - -void Solver::ResetSolution(){ - q.Reset(); - qdot.Reset(); - qddot.Reset(); - q_n.Reset(); - qdot_n.Reset(); - qddot_n.Reset(); -} - -void Solver::SaveToThePast(){ - q_n = q; - qdot_n = qdot; - qddot_n = qddot; - Loads_n = Loads; -} - -void Solver::SetInitialState(Config *config, Structure* structure){} - -/*CLASS ALPHAGENSOLVER*/ -AlphaGenSolver::AlphaGenSolver(unsigned int nDof, double val_rho, bool bool_linear) : Solver(nDof, bool_linear) { - - a.Initialize(nDof, 0.0); - a_n.Initialize(nDof, 0.0); - - rho = val_rho; - alpha_m = (2*rho-1)/(rho+1); - alpha_f = rho/(rho+1); - gamma = 0.5+alpha_f-alpha_m; - beta = 0.25*pow((gamma+0.5),2); - cout << "Integration with the alpha-generalized algorithm :" << endl; - cout << "rho : " << rho << endl; - cout << "alpha_m : " << alpha_m << endl; - cout << "alpha_f : " << alpha_f << endl; - cout << "gamma : " << gamma << endl; - cout << "beta : " << beta << endl; -} - -AlphaGenSolver::~AlphaGenSolver() {} - -CVector & AlphaGenSolver::GetAccVar(){ - return a; -} - -CVector & AlphaGenSolver::GetAccVar_n(){ - return a_n; -} - -void AlphaGenSolver::Iterate(double &t0, double &tf, Structure *structure){ - - double deltaT(tf-t0), epsilon(1e-6); - int nMaxIter(1000), nIter(0); - - gammaPrime = gamma/(deltaT*beta); - betaPrime = (1-alpha_m)/(pow(deltaT,2)*beta*(1-alpha_f)); - - /*--- Prediction phase ---*/ - qddot.Reset(); - a.Reset(); - - a += ScalVecProd(alpha_f/(1-alpha_m),qddot_n); - a -= ScalVecProd(alpha_m/(1-alpha_m),a_n); - - q = q_n; - q += ScalVecProd(deltaT,qdot_n); - q += ScalVecProd((0.5-beta)*deltaT*deltaT,a_n); - q += ScalVecProd(deltaT*deltaT*beta,a); - - qdot = qdot_n; - qdot += ScalVecProd((1-gamma)*deltaT,a_n); - qdot += ScalVecProd(deltaT*gamma,a); - - /*--- Tangent operator and corrector computation ---*/ - CVector res(qddot.GetSize(), 0.0); - CVector Deltaq(qddot.GetSize(), 0.0); - CMatrix St(qddot.GetSize(), qddot.GetSize(), 0.0); - ComputeResidual(structure,res); - while (res.norm() >= epsilon && nIter < nMaxIter){ - St.Reset(); - ComputeTangentOperator(structure,St); - SolveSys(St,res); - //*res -= ScalVecProd(double(2),res); //=deltaq - Deltaq.Reset(); - Deltaq += ScalVecProd(-1,res); - q += Deltaq; - qdot += ScalVecProd(gammaPrime,Deltaq); - qddot += ScalVecProd(betaPrime,Deltaq); - res.Reset(); - ComputeResidual(structure,res); - nIter++; - } - a += ScalVecProd((1-alpha_f)/(1-alpha_m),qddot); - -} - -void AlphaGenSolver::ComputeRHS(Structure* structure, CVector &RHS){ - - unsigned long size = q.GetSize(); - CMatrix CC(size, size, 0.0); - CMatrix KK(size, size, 0.0); - CVector NonLinTerm(size, 0.0); - - double cos_a; - - - if(structure->GetnDof() == 1){ - KK.SetElm(1,1,structure->Get_Kh()); - CC.SetElm(1,1,structure->Get_Ch()); - } - else if(structure->GetnDof() == 2){ - if(linear){ - cos_a = 1.0; - } - else{ - cos_a = cos(q[1]); - NonLinTerm[0] = (structure->Get_S())*sin(q[1])*pow(qdot[1],2); - } - CC.SetElm(1,1,structure->Get_Ch()); - CC.SetElm(2,2,structure->Get_Ca()); - KK.SetElm(1,1,structure->Get_Kh()); - KK.SetElm(2,2,structure->Get_Ka()); - } - else{ - - } - - RHS += Loads; - RHS -= MatVecProd(CC, qdot); - RHS -= MatVecProd(KK, q); - RHS += NonLinTerm; -} - -void AlphaGenSolver::ComputeResidual(Structure *structure, CVector & res){ - - res.Reset(); - - unsigned long size = q.GetSize(); - CMatrix MM(size, size, 0.0); - double cos_a; - - - if(structure->GetnDof() == 1){ - MM.SetElm(1,1, structure->Get_m()); - } - else if(structure->GetnDof() == 2){ - if(linear){ - cos_a = 1.0; - } - else{ - cos_a = cos(q[1]); - } - MM.SetElm(1,1, structure->Get_m()); - MM.SetElm(1,2,(structure->Get_S())*cos_a); - MM.SetElm(2,1,(structure->Get_S())*cos_a); - MM.SetElm(2,2,structure->Get_If()); - } - else{ - - } - - CVector RHS(size, 0.0); - ComputeRHS(structure, RHS); - - res = MatVecProd(MM, qddot) - RHS; -} - -void AlphaGenSolver::ComputeTangentOperator(Structure* structure, CMatrix &St){ - - St.Reset(); - - unsigned long size = q.GetSize(); - CMatrix MM(size, size, 0.0); - CMatrix Ct(size, size, 0.0); - CMatrix Kt(size, size, 0.0); - - if(structure->GetnDof() == 1){ - MM.SetElm(1,1, structure->Get_m()); - Kt.SetElm(1,1,(structure->Get_Kh())); - Ct.SetElm(1,1,(structure->Get_Ch())); - } - else if(structure->GetnDof() == 2){ - if(linear){ - MM.SetElm(1,1, structure->Get_m()); - MM.SetElm(1,2,(structure->Get_S())); - MM.SetElm(2,1,(structure->Get_S())); - MM.SetElm(2,2,structure->Get_If()); - Ct.SetElm(1,1,(structure->Get_Ch())); - Ct.SetElm(2,2,(structure->Get_Ca())); - Kt.SetElm(1,1,(structure->Get_Kh())); - Kt.SetElm(2,2,(structure->Get_Ka())); - } - else{ - MM.SetElm(1,1, structure->Get_m()); - MM.SetElm(1,2,(structure->Get_S())*cos(q[1])); - MM.SetElm(2,1,(structure->Get_S())*cos(q[1])); - MM.SetElm(2,2,structure->Get_If()); - Ct.SetElm(1,1,-(structure->Get_Ch())); - Ct.SetElm(1,2,(structure->Get_S())*sin(q[1])*2*qdot[1]); - Ct.SetElm(2,2,-(structure->Get_Ca())); - Kt.SetElm(1,1,-(structure->Get_Kh())); - Kt.SetElm(1,2,(structure->Get_S())*cos(q[1])*pow(qdot[1],2)); - Kt.SetElm(2,2,-(structure->Get_Ka())); - } - - } - else{ - - } - - St += ScalMatProd(betaPrime, MM); - St += ScalMatProd(gammaPrime, Ct); - St += Kt; - -} - -void AlphaGenSolver::ResetSolution(){ - - Solver::ResetSolution(); - a.Reset(); - a_n.Reset(); -} - -void AlphaGenSolver::SaveToThePast(){ - - Solver::SaveToThePast(); - a_n = a; - -} - -void AlphaGenSolver::SetInitialState(Config* config, Structure* structure){ - - if(config->GetRestartSol() == "YES"){ - string InputFileName = config->GetRestartFile(); - string text_line; - string token, tempString; - size_t pos; - string delimiter = "\t"; - ifstream InputFile; - InputFile.open(InputFileName.c_str(), ios::in); - double buffer[(4*structure->GetnDof())+1]; - int kk = 0; - int jj; - while (getline(InputFile,text_line)){ - tempString = text_line; - jj = 0; - if (kk == 1){ - while ((pos = tempString.find(delimiter)) != string::npos){ - token = tempString.substr(0,pos); - tempString.erase(0,pos+delimiter.length()); - buffer[jj] = atof(token.c_str()); - jj += 1; - } - buffer[jj] = atof(tempString.c_str()); - - if(structure->GetnDof() == 1){ - q_n[0] = buffer[1]; - qdot_n[0] = buffer[2]; - qddot_n[0] = buffer[3]; - a_n[0] = buffer[4]; - } - else if (structure->GetnDof() == 2){ - q_n[0] = buffer[1]; - q_n[1] = buffer[2]; - qdot_n[0] = buffer[3]; - qdot_n[1] = buffer[4]; - qddot_n[0] = buffer[5]; - qddot_n[1] = buffer[6]; - a_n[0] = buffer[7]; - a_n[1] = buffer[8]; - } - q_n.print(); - qdot_n.print(); - qddot_n.print(); - a_n.print(); - } - else if (kk == 2){ - while ((pos = tempString.find(delimiter)) != string::npos){ - token = tempString.substr(0,pos); - tempString.erase(0,pos+delimiter.length()); - buffer[jj] = atof(token.c_str()); - jj += 1; - } - buffer[jj] = atof(tempString.c_str()); - - if(structure->GetnDof() == 1){ - q[0] = buffer[1]; - qdot[0] = buffer[2]; - qddot[0] = buffer[3]; - a[0] = buffer[4]; - } - else if (structure->GetnDof() == 2){ - q[0] = buffer[1]; - q[1] = buffer[2]; - qdot[0] = buffer[3]; - qdot[1] = buffer[4]; - qddot[0] = buffer[5]; - qddot[1] = buffer[6]; - a[0] = buffer[7]; - a[1] = buffer[8]; - } - q.print(); - qdot.print(); - qddot.print(); - a.print(); - } - kk += 1; - } - InputFile.close(); - } - else{ - cout << "Setting basic initial conditions for alpha-Gen" << endl; - q.Reset(); - q_n.Reset(); - cout << "Read initial configuration" << endl; - q[0] = config->GetInitialDisp(); - if(structure->GetnDof() == 2) q[1] = config->GetInitialAngle(); - cout << "Initial plunging displacement : " << q[0] << endl; - cout << "Initial pitching displacement : " << q[1] << endl; - - qdot.Reset(); - qddot.Reset(); - - unsigned long size = q.GetSize(); - CVector RHS(size, 0.0); - CMatrix MM(size, size, 0.0); - - double cos_a; - - - if(structure->GetnDof() == 1){ - MM.SetElm(1,1, structure->Get_m()); - } - else if(structure->GetnDof() == 2){ - if(linear){ - cos_a = 1.0; - } - else{ - cos_a = cos(q[1]); - } - MM.SetElm(1,1, structure->Get_m()); - MM.SetElm(1,2,(structure->Get_S())*cos_a); - MM.SetElm(2,1,(structure->Get_S())*cos_a); - MM.SetElm(2,2,structure->Get_If()); - } - else{ - - } - - ComputeRHS(structure, RHS); - SolveSys(MM, RHS); - qddot = RHS; - a = qddot; - - } -} - -/*CLASS RK4 SOLVER*/ -RK4Solver::RK4Solver(unsigned nDof, bool bool_linear) : Solver(nDof, bool_linear){ - size = nDof; - lastTime = 0.0; - currentTime = 0.0; -} - -RK4Solver::~RK4Solver(){} - -void RK4Solver::Iterate(double& t0, double& tf, Structure *structure){ - - double h = tf-t0; - lastTime = t0; - currentTime = tf; - - CVector k1(2*size); - CVector k2(2*size); - CVector k3(2*size); - CVector k4(2*size); - - CVector state0(2*size); - CVector statef(2*size); - CVector statef_dot(2*size); - state0 = SetState_n(); - - CVector TEMP(2*size); - - EvaluateStateDerivative(lastTime, state0, k1, structure); - TEMP = state0+(k1*(h/2.0)); - EvaluateStateDerivative(lastTime+h/2.0, TEMP, k2, structure); - TEMP = state0+(k2*(h/2.0)); - EvaluateStateDerivative(lastTime+h/2.0, TEMP, k3, structure); - TEMP = state0+(k3*h); - EvaluateStateDerivative(lastTime+h, TEMP, k4, structure); - - statef = state0 + ((k1 + k2*2.0 + k3*2.0 + k4)*(h/6.0)); - - EvaluateStateDerivative(lastTime+h, statef, statef_dot, structure); - - if(structure->GetnDof() == 1){ - q[0] = statef[0]; - qdot[0] = statef[1]; - qddot[0] = statef_dot[1]; - } - else if (structure->GetnDof() == 2){ - q[0] = statef[0]; - q[1] = statef[1]; - qdot[0] = statef[2]; - qdot[1] = statef[3]; - qddot[0] = statef_dot[2]; - qddot[1] = statef_dot[3]; - } - else{ - - } -} - -void RK4Solver::EvaluateStateDerivative(double tCurrent, CVector &state, CVector &stateDerivative, Structure* structure){ - - CVector stateLoads(size, 0.0); - interpLoads(tCurrent, stateLoads); - - CMatrix MM(size, size, 0.0); - CMatrix CC(size, size, 0.0); - CMatrix KK(size, size, 0.0); - CVector NonLinTerm(size, 0.0); - CVector RHS(size, 0.0); - - CVector q_current(size, 0.0); - CVector qdot_current(size, 0.0); - CVector qddot_current(size, 0.0); - - if(structure->GetnDof() == 1){ - MM.SetElm(1,1, structure->Get_m()); - KK.SetElm(1,1,structure->Get_Kh()); - CC.SetElm(1,1,structure->Get_Ch()); - q_current[0] = state[0]; - qdot_current[0] = state[1]; - stateDerivative[0] = state[1]; - } - else if (structure->GetnDof() == 2){ - double cos_a; - if(linear){ - cos_a = 1.0; - } - else{ - cos_a = cos(state[1]); - NonLinTerm[0] = (structure->Get_S())*sin(state[1])*pow(state[3],2); - } - MM.SetElm(1,1, structure->Get_m()); - MM.SetElm(1,2,(structure->Get_S())*cos_a); - MM.SetElm(2,1,(structure->Get_S())*cos_a); - MM.SetElm(2,2,structure->Get_If()); - CC.SetElm(1,1,structure->Get_Ch()); - CC.SetElm(2,2,structure->Get_Ca()); - KK.SetElm(1,1,structure->Get_Kh()); - KK.SetElm(2,2,structure->Get_Ka()); - q_current[0] = state[0]; - q_current[1] = state[1]; - qdot_current[0] = state[2]; - qdot_current[1] = state[3]; - stateDerivative[0] = state[2]; - stateDerivative[1] = state[3]; - } - else{ - - } - - RHS += stateLoads; - RHS -= MatVecProd(CC, qdot_current); - RHS -= MatVecProd(KK, q_current); - RHS += NonLinTerm; - - SolveSys(MM, RHS); - qddot_current = RHS; - - if(structure->GetnDof() == 1){ - stateDerivative[1] = qddot_current[0]; - } - else if (structure->GetnDof() == 2){ - stateDerivative[2] = qddot_current[0]; - stateDerivative[3] = qddot_current[1]; - } - else{ - - } - -} - -void RK4Solver::interpLoads(double &tCurrent, CVector &val_loads){ - - if (lastTime != currentTime){ - val_loads[0] = (Loads[0] - Loads_n[0])/(currentTime-lastTime)*(tCurrent - lastTime) + Loads_n[0]; - if(size == 2) val_loads[1] = (Loads[1] - Loads_n[1])/(currentTime-lastTime)*(tCurrent - lastTime) + Loads_n[1]; - } - else{ - val_loads[0] = Loads[0]; - if(size == 2) val_loads[1] = Loads[1]; - } -} - -void RK4Solver::SetInitialState(Config *config, Structure *structure){ - - if (config->GetRestartSol() == "YES"){ - - } - else{ - cout << "Setting basic initial conditions for RK4" << endl; - q.Reset(); - q_n.Reset(); - cout << "Read initial configuration" << endl; - q[0] = config->GetInitialDisp(); - if(structure->GetnDof() == 2) q[1] = config->GetInitialAngle(); - cout << "Initial plunging displacement : " << q[0] << endl; - cout << "Initial pitching displacement : " << q[1] << endl; - qdot.Reset(); - - lastTime = 0.0; - currentTime = 0.0; - - qddot.Reset(); - CVector state(2*size); - CVector state_dot(2*size); - - state = SetState(); - EvaluateStateDerivative(0.0, state, state_dot, structure); - - if(structure->GetnDof() == 1){ - qddot[0] = state_dot[1]; - } - else if (structure->GetnDof() == 2){ - qddot[0] = state_dot[2]; - qddot[1] = state_dot[3]; - } - else{ - - } - - } -} - -CVector RK4Solver::SetState(){ - - CVector state(2*size); - - if(size == 1){ - state[0] = q[0]; - state[1] = qdot[0]; - } - else if (size == 2){ - state[0] = q[0]; - state[1] = q[1]; - state[2] = qdot[0]; - state[3] = qdot[1]; - } - else{ - - } - - return state; - -} - -CVector RK4Solver::SetState_n(){ - - CVector state(2*size); - - if(size == 1){ - state[0] = q_n[0]; - state[1] = qdot_n[0]; - } - else if (size == 2){ - state[0] = q_n[0]; - state[1] = q_n[1]; - state[2] = qdot_n[0]; - state[3] = qdot_n[1]; - } - else{ - - } - - return state; -} diff --git a/src/structure.cpp b/src/structure.cpp deleted file mode 100644 index fedd2e5..0000000 --- a/src/structure.cpp +++ /dev/null @@ -1,172 +0,0 @@ -#include "../include/structure.h" -#include -#include -#include -#include -#include - - -using namespace std; - -Structure::Structure(Config* config){ - - centerOfRotation[0] = 0.0; - centerOfRotation[1] = 0.0; - centerOfRotation[2] = 0.0; - - centerOfRotation_n[0] = 0.0; - centerOfRotation_n[1] = 0.0; - centerOfRotation_n[2] = 0.0; - - - if(config->GetStructType() == "SPRING_HOR" || config->GetStructType() == "SPRING_VER" ){ - nDof = 1; - m = config->GetSpringMass(); - Kh = config->GetSpringStiffness(); - Ch = config->GetSpringDamping(); - Ka = 0; - Ca = 0; - xf = 0; - xCG = 0; - c = 0; - b = 0; - S = 0; - ICG = 0; - If = 0; - cout << "Setting mass-spring-damper system" << endl; - cout << "Number of DOF : " << nDof << endl; - cout << "Plunging mass : " << m << " [kg]" << endl; - cout << "Plunging damping : " << Ch << " [Ns/m]" << endl; - cout << "Plunging stiffness : " << Kh << " [N/m]" << endl; - } - else if (config->GetStructType() == "AIRFOIL"){ - nDof = 2; - m = config->GetSpringMass(); - Kh = config->GetSpringStiffness(); - Ch = config->GetSpringDamping(); - Ka = config->GetTorsionalStiffness(); - Ca = config->GetTorsionalDamping(); - xf = config->GetFlexuralAxis(); - xCG = config->GetGravityCenter(); - c = config->GetCord(); - b = c/2.0; - S = m*(xCG-xf); - //If = ICG + m*pow((xCG-xf),2); - If = config->GetInertiaFlexural(); - //S = m*(c/2.0-xf); - //Ia = 1.0/3.0*m*(c*c-3*c*xf+3*xf*xf); - cout << "Setting pitching-plunging airfoil system" << endl; - cout << "Number of DOF : " << nDof << endl; - cout << "Airfoil mass : " << m << " [kg]" << endl; - cout << "Airfoil cord : " << c << " [m]" << endl; - cout << "Position of the flexural axis : " << xf << " [m]" << endl; - cout << "Inertia around the flexural axis : " << If << " [kg m²]" << endl; - cout << "Plunging damping : " << Ch << " [Ns/m]" << endl; - cout << "Plunging stiffness : " << Kh << " [N/m]" << endl; - cout << "Pitching damping : " << Ca << " [Ns]" << endl; - cout << "Pitching stiffness : " << Ka << " [N]" << endl; - cout << "Position of the center of gravity : " << xCG << " [m]" << endl; - cout << "Static unbalance : " << S << " [kg m]" << endl; - - centerOfRotation[0] = xf; - centerOfRotation_n[0] = xf; - } - else{ - nDof = 0; - cerr << "Invalid structural type. Available choices are : SPRIN_HOR, SPRING_VER and AIRFOIL." << endl; - throw(-1); - } -} - -Structure::~Structure(){} - -void Structure::SetCenterOfRotation_X(double coord_x){ - centerOfRotation[0] = coord_x; -} - -void Structure::SetCenterOfRotation_Y(double coord_y){ - centerOfRotation[1] = coord_y; -} - -void Structure::SetCenterOfRotation_Z(double coord_z){ - centerOfRotation[2] = coord_z; -} - -double Structure::GetCenterOfRotation_x() const{ - return centerOfRotation[0]; -} - -double Structure::GetCenterOfRotation_y() const{ - return centerOfRotation[1]; -} - -double Structure::GetCenterOfRotation_z() const{ - return centerOfRotation[2]; -} - -const double* Structure::GetCenterOfRotation() const{ - return centerOfRotation; -} - -void Structure::SetCenterOfRotation_n_X(double coord_x){ - centerOfRotation_n[0] = coord_x; -} - -void Structure::SetCenterOfRotation_n_Y(double coord_y){ - centerOfRotation_n[1] = coord_y; -} - -void Structure::SetCenterOfRotation_n_Z(double coord_z){ - centerOfRotation_n[2] = coord_z; -} - -double Structure::GetCenterOfRotation_n_x() const{ - return centerOfRotation_n[0]; -} - -double Structure::GetCenterOfRotation_n_y() const{ - return centerOfRotation_n[1]; -} - -double Structure::GetCenterOfRotation_n_z() const{ - return centerOfRotation_n[2]; -} - -unsigned int Structure::GetnDof() const{ - return nDof; -} - -double Structure::Get_m() const{ - - return m; -} - -double Structure::Get_Kh() const{ - - return Kh; -} - -double Structure::Get_Ka() const{ - - return Ka; -} - -double Structure::Get_Ch() const{ - - return Ch; -} - -double Structure::Get_Ca() const{ - - return Ca; -} - -double Structure::Get_S() const{ - - return S; -} - -double Structure::Get_If() const{ - - return If; -} diff --git a/src/testcvector.cpp b/src/testcvector.cpp deleted file mode 100644 index 9773ec5..0000000 --- a/src/testcvector.cpp +++ /dev/null @@ -1,122 +0,0 @@ -#include - -#include "../include/MatVec.h" - -int main(int argc, char** argv){ - - //Test initialization - CVector myVec(5, 5.0); - - //Test copy - CVector myVec2(myVec); - - //Test initiate - std::cout << "initiate()"<< std::endl; - CVector myVec4; - myVec4.Initialize(5,7.0); - myVec4.print(); - - //Test getSize() - std::cout << "getSize()"<< std::endl; - std::cout << myVec.GetSize() << std::endl; - std::cout << myVec2.GetSize() << std::endl; - - //Test dotProduct() - std::cout << "dot()"<< std::endl; - CVector myVec3(5,1.0); - double dotProduct(0.0), dotProduct2(0.0); - dotProduct = myVec.dotProd(myVec3); - dotProduct2 = myVec3.dotProd(myVec); - std::cout << dotProduct << std::endl; - std::cout << dotProduct2 << std::endl; - - //Test norm() - std::cout << "norm()"<< std::endl; - std::cout << myVec.norm() << std::endl; - - //Test display() - std::cout << "display()"<< std::endl; - myVec.print(); - - //Test reset(); - std::cout << "reset()"<< std::endl; - myVec2.Reset(); - myVec2.print(); - - //Test operator[] as member - std::cout << "operator[]"<< std::endl; - std::cout << myVec[3] << std::endl; - myVec[3] = 100.0; - myVec.print(); - - //Test operator= as member - std::cout << "operator="<< std::endl; - myVec2 = myVec; - myVec2.print(); - - //Test operator += as member - std::cout << "operator+="<< std::endl; - myVec.SetAllValues(2.0); - myVec2.SetAllValues(3.0); - myVec2 += myVec; - myVec2.print(); - - //Test operator -= as member - std::cout << "operator-="<< std::endl; - myVec2 -= myVec; - myVec2.print(); - - //Test operator *= as member - std::cout << "operator*="<< std::endl; - myVec.SetAllValues(2.0); - myVec *= 2.0; - myVec.print(); - - //Test operator /= as member - std::cout << "operator/="<< std::endl; - myVec /= 2.0; - myVec.print(); - - //Test addition - std::cout << "addition"<< std::endl; - myVec.SetAllValues(2.0); - myVec2.SetAllValues(3.0); - CVector vecSum(5, 0.0); - vecSum = myVec + myVec2; - vecSum.print(); - - //Test multiplication - std::cout << "multiplication"<< std::endl; - myVec.SetAllValues(2.0); - CVector vecMul(5, 0.0); - vecMul = myVec*2.0; - vecMul = 3.0*myVec; - vecMul.print(); - - //Test linear combination - std::cout << "linear combination"<< std::endl; - myVec.SetAllValues(1.0); - myVec2.SetAllValues(1.0); - myVec3.SetAllValues(1.0); - myVec4 = myVec + (2.0*myVec2 + 5.0*myVec3)*2.0; - myVec4.print(); - - std::cout << "Pointers..."<< std::endl; - CVector* myVecPoint = NULL; - CVector* myVecPoint2 = NULL; - myVecPoint = new CVector(5,10.0); - myVecPoint->print(); - std::cout << (*myVecPoint)[3] << std::endl; - myVecPoint2 = new CVector(*myVecPoint); - myVecPoint2->print(); - myVecPoint->Reset(); - *myVecPoint2 = *myVecPoint; - myVecPoint2->print(); - - - if(myVecPoint != NULL) delete myVecPoint; - if(myVecPoint2 != NULL) delete myVecPoint2; - - return 0; -} -