diff options
38 files changed, 790 insertions, 3570 deletions
diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000..6080600 --- /dev/null +++ b/.gitmodules @@ -0,0 +1,3 @@ +[submodule "src/randutils"] + path = src/randutils + url = https://gist.github.com/imneme/540829265469e673d045 diff --git a/CMakeLists.txt b/CMakeLists.txt index e924e98..6746aa9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,20 +1,14 @@ -cmake_minimum_required(VERSION 3.5) -project(fracture) +cmake_minimum_required(VERSION 3.9) -include_directories(src ~/.local/include /usr/include/suitesparse) -link_directories(~/.local/lib) +set(CMAKE_CXX_FLAGS_DEBUG "-g -Wall") +set(CMAKE_CXX_FLAGS_RELEASE "-O3 -flto") -file(GLOB LIB_SOURCES lib/*.c) +set (CMAKE_CXX_STANDARD 17) +set (CMAKE_C_STANDARD 11) -add_library(fracture_stuff ${LIB_SOURCES}) +include(GNUInstallDirs) -file(GLOB EXE_SOURCES src/*.c) - -foreach( src_file ${EXE_SOURCES} ) - string( REGEX REPLACE ".*/src/(.*)\.c" "\\1" exe_name ${src_file} ) - add_executable( ${exe_name} ${src_file} ) - target_link_libraries(${exe_name} fracture_stuff gsl c cblas lapack dl pthread cholmod amd colamd suitesparseconfig camd ccolamd rt metis fftw3 m jst) - install(TARGETS ${exe_name} DESTINATION bin) -endforeach( src_file ${EXE_SOURCES} ) +add_subdirectory(lib) +add_subdirectory(src) diff --git a/LICENSE b/LICENSE deleted file mode 100644 index 0d1fadb..0000000 --- a/LICENSE +++ /dev/null @@ -1,674 +0,0 @@ - GNU GENERAL PUBLIC LICENSE - Version 3, 29 June 2007 - - Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/> - Everyone is permitted to copy and distribute verbatim copies - of this license document, but changing it is not allowed. - - Preamble - - The GNU General Public License is a free, copyleft license for -software and other kinds of works. - - The licenses for most software and other practical works are designed -to take away your freedom to share and change the works. By contrast, -the GNU General Public License is intended to guarantee your freedom to -share and change all versions of a program--to make sure it remains free -software for all its users. We, the Free Software Foundation, use the -GNU General Public License for most of our software; it applies also to -any other work released this way by its authors. You can apply it to -your programs, too. - - When we speak of free software, we are referring to freedom, not -price. Our General Public Licenses are designed to make sure that you -have the freedom to distribute copies of free software (and charge for -them if you wish), that you receive source code or can get it if you -want it, that you can change the software or use pieces of it in new -free programs, and that you know you can do these things. - - To protect your rights, we need to prevent others from denying you -these rights or asking you to surrender the rights. Therefore, you have -certain responsibilities if you distribute copies of the software, or if -you modify it: responsibilities to respect the freedom of others. - - For example, if you distribute copies of such a program, whether -gratis or for a fee, you must pass on to the recipients the same -freedoms that you received. You must make sure that they, too, receive -or can get the source code. And you must show them these terms so they -know their rights. - - Developers that use the GNU GPL protect your rights with two steps: -(1) assert copyright on the software, and (2) offer you this License -giving you legal permission to copy, distribute and/or modify it. - - For the developers' and authors' protection, the GPL clearly explains -that there is no warranty for this free software. For both users' and -authors' sake, the GPL requires that modified versions be marked as -changed, so that their problems will not be attributed erroneously to -authors of previous versions. - - Some devices are designed to deny users access to install or run -modified versions of the software inside them, although the manufacturer -can do so. This is fundamentally incompatible with the aim of -protecting users' freedom to change the software. The systematic -pattern of such abuse occurs in the area of products for individuals to -use, which is precisely where it is most unacceptable. Therefore, we -have designed this version of the GPL to prohibit the practice for those -products. If such problems arise substantially in other domains, we -stand ready to extend this provision to those domains in future versions -of the GPL, as needed to protect the freedom of users. - - Finally, every program is threatened constantly by software patents. -States should not allow patents to restrict development and use of -software on general-purpose computers, but in those that do, we wish to -avoid the special danger that patents applied to a free program could -make it effectively proprietary. To prevent this, the GPL assures that -patents cannot be used to render the program non-free. - - The precise terms and conditions for copying, distribution and -modification follow. - - TERMS AND CONDITIONS - - 0. Definitions. - - "This License" refers to version 3 of the GNU General Public License. - - "Copyright" also means copyright-like laws that apply to other kinds of -works, such as semiconductor masks. - - "The Program" refers to any copyrightable work licensed under this -License. Each licensee is addressed as "you". "Licensees" and -"recipients" may be individuals or organizations. - - To "modify" a work means to copy from or adapt all or part of the work -in a fashion requiring copyright permission, other than the making of an -exact copy. The resulting work is called a "modified version" of the -earlier work or a work "based on" the earlier work. - - A "covered work" means either the unmodified Program or a work based -on the Program. - - To "propagate" a work means to do anything with it that, without -permission, would make you directly or secondarily liable for -infringement under applicable copyright law, except executing it on a -computer or modifying a private copy. Propagation includes copying, -distribution (with or without modification), making available to the -public, and in some countries other activities as well. - - To "convey" a work means any kind of propagation that enables other -parties to make or receive copies. Mere interaction with a user through -a computer network, with no transfer of a copy, is not conveying. - - An interactive user interface displays "Appropriate Legal Notices" -to the extent that it includes a convenient and prominently visible -feature that (1) displays an appropriate copyright notice, and (2) -tells the user that there is no warranty for the work (except to the -extent that warranties are provided), that licensees may convey the -work under this License, and how to view a copy of this License. If -the interface presents a list of user commands or options, such as a -menu, a prominent item in the list meets this criterion. - - 1. Source Code. - - The "source code" for a work means the preferred form of the work -for making modifications to it. "Object code" means any non-source -form of a work. - - A "Standard Interface" means an interface that either is an official -standard defined by a recognized standards body, or, in the case of -interfaces specified for a particular programming language, one that -is widely used among developers working in that language. - - The "System Libraries" of an executable work include anything, other -than the work as a whole, that (a) is included in the normal form of -packaging a Major Component, but which is not part of that Major -Component, and (b) serves only to enable use of the work with that -Major Component, or to implement a Standard Interface for which an -implementation is available to the public in source code form. A -"Major Component", in this context, means a major essential component -(kernel, window system, and so on) of the specific operating system -(if any) on which the executable work runs, or a compiler used to -produce the work, or an object code interpreter used to run it. - - The "Corresponding Source" for a work in object code form means all -the source code needed to generate, install, and (for an executable -work) run the object code and to modify the work, including scripts to -control those activities. However, it does not include the work's -System Libraries, or general-purpose tools or generally available free -programs which are used unmodified in performing those activities but -which are not part of the work. For example, Corresponding Source -includes interface definition files associated with source files for -the work, and the source code for shared libraries and dynamically -linked subprograms that the work is specifically designed to require, -such as by intimate data communication or control flow between those -subprograms and other parts of the work. - - The Corresponding Source need not include anything that users -can regenerate automatically from other parts of the Corresponding -Source. - - The Corresponding Source for a work in source code form is that -same work. - - 2. Basic Permissions. - - All rights granted under this License are granted for the term of -copyright on the Program, and are irrevocable provided the stated -conditions are met. This License explicitly affirms your unlimited -permission to run the unmodified Program. The output from running a -covered work is covered by this License only if the output, given its -content, constitutes a covered work. This License acknowledges your -rights of fair use or other equivalent, as provided by copyright law. - - You may make, run and propagate covered works that you do not -convey, without conditions so long as your license otherwise remains -in force. You may convey covered works to others for the sole purpose -of having them make modifications exclusively for you, or provide you -with facilities for running those works, provided that you comply with -the terms of this License in conveying all material for which you do -not control copyright. Those thus making or running the covered works -for you must do so exclusively on your behalf, under your direction -and control, on terms that prohibit them from making any copies of -your copyrighted material outside their relationship with you. - - Conveying under any other circumstances is permitted solely under -the conditions stated below. Sublicensing is not allowed; section 10 -makes it unnecessary. - - 3. Protecting Users' Legal Rights From Anti-Circumvention Law. - - No covered work shall be deemed part of an effective technological -measure under any applicable law fulfilling obligations under article -11 of the WIPO copyright treaty adopted on 20 December 1996, or -similar laws prohibiting or restricting circumvention of such -measures. - - When you convey a covered work, you waive any legal power to forbid -circumvention of technological measures to the extent such circumvention -is effected by exercising rights under this License with respect to -the covered work, and you disclaim any intention to limit operation or -modification of the work as a means of enforcing, against the work's -users, your or third parties' legal rights to forbid circumvention of -technological measures. - - 4. Conveying Verbatim Copies. - - You may convey verbatim copies of the Program's source code as you -receive it, in any medium, provided that you conspicuously and -appropriately publish on each copy an appropriate copyright notice; -keep intact all notices stating that this License and any -non-permissive terms added in accord with section 7 apply to the code; -keep intact all notices of the absence of any warranty; and give all -recipients a copy of this License along with the Program. - - You may charge any price or no price for each copy that you convey, -and you may offer support or warranty protection for a fee. - - 5. Conveying Modified Source Versions. - - You may convey a work based on the Program, or the modifications to -produce it from the Program, in the form of source code under the -terms of section 4, provided that you also meet all of these conditions: - - a) The work must carry prominent notices stating that you modified - it, and giving a relevant date. - - b) The work must carry prominent notices stating that it is - released under this License and any conditions added under section - 7. This requirement modifies the requirement in section 4 to - "keep intact all notices". - - c) You must license the entire work, as a whole, under this - License to anyone who comes into possession of a copy. This - License will therefore apply, along with any applicable section 7 - additional terms, to the whole of the work, and all its parts, - regardless of how they are packaged. This License gives no - permission to license the work in any other way, but it does not - invalidate such permission if you have separately received it. - - d) If the work has interactive user interfaces, each must display - Appropriate Legal Notices; however, if the Program has interactive - interfaces that do not display Appropriate Legal Notices, your - work need not make them do so. - - A compilation of a covered work with other separate and independent -works, which are not by their nature extensions of the covered work, -and which are not combined with it such as to form a larger program, -in or on a volume of a storage or distribution medium, is called an -"aggregate" if the compilation and its resulting copyright are not -used to limit the access or legal rights of the compilation's users -beyond what the individual works permit. Inclusion of a covered work -in an aggregate does not cause this License to apply to the other -parts of the aggregate. - - 6. Conveying Non-Source Forms. - - You may convey a covered work in object code form under the terms -of sections 4 and 5, provided that you also convey the -machine-readable Corresponding Source under the terms of this License, -in one of these ways: - - a) Convey the object code in, or embodied in, a physical product - (including a physical distribution medium), accompanied by the - Corresponding Source fixed on a durable physical medium - customarily used for software interchange. - - b) Convey the object code in, or embodied in, a physical product - (including a physical distribution medium), accompanied by a - written offer, valid for at least three years and valid for as - long as you offer spare parts or customer support for that product - model, to give anyone who possesses the object code either (1) a - copy of the Corresponding Source for all the software in the - product that is covered by this License, on a durable physical - medium customarily used for software interchange, for a price no - more than your reasonable cost of physically performing this - conveying of source, or (2) access to copy the - Corresponding Source from a network server at no charge. - - c) Convey individual copies of the object code with a copy of the - written offer to provide the Corresponding Source. This - alternative is allowed only occasionally and noncommercially, and - only if you received the object code with such an offer, in accord - with subsection 6b. - - d) Convey the object code by offering access from a designated - place (gratis or for a charge), and offer equivalent access to the - Corresponding Source in the same way through the same place at no - further charge. You need not require recipients to copy the - Corresponding Source along with the object code. If the place to - copy the object code is a network server, the Corresponding Source - may be on a different server (operated by you or a third party) - that supports equivalent copying facilities, provided you maintain - clear directions next to the object code saying where to find the - Corresponding Source. Regardless of what server hosts the - Corresponding Source, you remain obligated to ensure that it is - available for as long as needed to satisfy these requirements. - - e) Convey the object code using peer-to-peer transmission, provided - you inform other peers where the object code and Corresponding - Source of the work are being offered to the general public at no - charge under subsection 6d. - - A separable portion of the object code, whose source code is excluded -from the Corresponding Source as a System Library, need not be -included in conveying the object code work. - - A "User Product" is either (1) a "consumer product", which means any -tangible personal property which is normally used for personal, family, -or household purposes, or (2) anything designed or sold for incorporation -into a dwelling. In determining whether a product is a consumer product, -doubtful cases shall be resolved in favor of coverage. For a particular -product received by a particular user, "normally used" refers to a -typical or common use of that class of product, regardless of the status -of the particular user or of the way in which the particular user -actually uses, or expects or is expected to use, the product. A product -is a consumer product regardless of whether the product has substantial -commercial, industrial or non-consumer uses, unless such uses represent -the only significant mode of use of the product. - - "Installation Information" for a User Product means any methods, -procedures, authorization keys, or other information required to install -and execute modified versions of a covered work in that User Product from -a modified version of its Corresponding Source. The information must -suffice to ensure that the continued functioning of the modified object -code is in no case prevented or interfered with solely because -modification has been made. - - If you convey an object code work under this section in, or with, or -specifically for use in, a User Product, and the conveying occurs as -part of a transaction in which the right of possession and use of the -User Product is transferred to the recipient in perpetuity or for a -fixed term (regardless of how the transaction is characterized), the -Corresponding Source conveyed under this section must be accompanied -by the Installation Information. But this requirement does not apply -if neither you nor any third party retains the ability to install -modified object code on the User Product (for example, the work has -been installed in ROM). - - The requirement to provide Installation Information does not include a -requirement to continue to provide support service, warranty, or updates -for a work that has been modified or installed by the recipient, or for -the User Product in which it has been modified or installed. Access to a -network may be denied when the modification itself materially and -adversely affects the operation of the network or violates the rules and -protocols for communication across the network. - - Corresponding Source conveyed, and Installation Information provided, -in accord with this section must be in a format that is publicly -documented (and with an implementation available to the public in -source code form), and must require no special password or key for -unpacking, reading or copying. - - 7. Additional Terms. - - "Additional permissions" are terms that supplement the terms of this -License by making exceptions from one or more of its conditions. -Additional permissions that are applicable to the entire Program shall -be treated as though they were included in this License, to the extent -that they are valid under applicable law. If additional permissions -apply only to part of the Program, that part may be used separately -under those permissions, but the entire Program remains governed by -this License without regard to the additional permissions. - - When you convey a copy of a covered work, you may at your option -remove any additional permissions from that copy, or from any part of -it. (Additional permissions may be written to require their own -removal in certain cases when you modify the work.) You may place -additional permissions on material, added by you to a covered work, -for which you have or can give appropriate copyright permission. - - Notwithstanding any other provision of this License, for material you -add to a covered work, you may (if authorized by the copyright holders of -that material) supplement the terms of this License with terms: - - a) Disclaiming warranty or limiting liability differently from the - terms of sections 15 and 16 of this License; or - - b) Requiring preservation of specified reasonable legal notices or - author attributions in that material or in the Appropriate Legal - Notices displayed by works containing it; or - - c) Prohibiting misrepresentation of the origin of that material, or - requiring that modified versions of such material be marked in - reasonable ways as different from the original version; or - - d) Limiting the use for publicity purposes of names of licensors or - authors of the material; or - - e) Declining to grant rights under trademark law for use of some - trade names, trademarks, or service marks; or - - f) Requiring indemnification of licensors and authors of that - material by anyone who conveys the material (or modified versions of - it) with contractual assumptions of liability to the recipient, for - any liability that these contractual assumptions directly impose on - those licensors and authors. - - All other non-permissive additional terms are considered "further -restrictions" within the meaning of section 10. If the Program as you -received it, or any part of it, contains a notice stating that it is -governed by this License along with a term that is a further -restriction, you may remove that term. If a license document contains -a further restriction but permits relicensing or conveying under this -License, you may add to a covered work material governed by the terms -of that license document, provided that the further restriction does -not survive such relicensing or conveying. - - If you add terms to a covered work in accord with this section, you -must place, in the relevant source files, a statement of the -additional terms that apply to those files, or a notice indicating -where to find the applicable terms. - - Additional terms, permissive or non-permissive, may be stated in the -form of a separately written license, or stated as exceptions; -the above requirements apply either way. - - 8. Termination. - - You may not propagate or modify a covered work except as expressly -provided under this License. Any attempt otherwise to propagate or -modify it is void, and will automatically terminate your rights under -this License (including any patent licenses granted under the third -paragraph of section 11). - - However, if you cease all violation of this License, then your -license from a particular copyright holder is reinstated (a) -provisionally, unless and until the copyright holder explicitly and -finally terminates your license, and (b) permanently, if the copyright -holder fails to notify you of the violation by some reasonable means -prior to 60 days after the cessation. - - Moreover, your license from a particular copyright holder is -reinstated permanently if the copyright holder notifies you of the -violation by some reasonable means, this is the first time you have -received notice of violation of this License (for any work) from that -copyright holder, and you cure the violation prior to 30 days after -your receipt of the notice. - - Termination of your rights under this section does not terminate the -licenses of parties who have received copies or rights from you under -this License. If your rights have been terminated and not permanently -reinstated, you do not qualify to receive new licenses for the same -material under section 10. - - 9. Acceptance Not Required for Having Copies. - - You are not required to accept this License in order to receive or -run a copy of the Program. Ancillary propagation of a covered work -occurring solely as a consequence of using peer-to-peer transmission -to receive a copy likewise does not require acceptance. However, -nothing other than this License grants you permission to propagate or -modify any covered work. These actions infringe copyright if you do -not accept this License. Therefore, by modifying or propagating a -covered work, you indicate your acceptance of this License to do so. - - 10. Automatic Licensing of Downstream Recipients. - - Each time you convey a covered work, the recipient automatically -receives a license from the original licensors, to run, modify and -propagate that work, subject to this License. You are not responsible -for enforcing compliance by third parties with this License. - - An "entity transaction" is a transaction transferring control of an -organization, or substantially all assets of one, or subdividing an -organization, or merging organizations. If propagation of a covered -work results from an entity transaction, each party to that -transaction who receives a copy of the work also receives whatever -licenses to the work the party's predecessor in interest had or could -give under the previous paragraph, plus a right to possession of the -Corresponding Source of the work from the predecessor in interest, if -the predecessor has it or can get it with reasonable efforts. - - You may not impose any further restrictions on the exercise of the -rights granted or affirmed under this License. For example, you may -not impose a license fee, royalty, or other charge for exercise of -rights granted under this License, and you may not initiate litigation -(including a cross-claim or counterclaim in a lawsuit) alleging that -any patent claim is infringed by making, using, selling, offering for -sale, or importing the Program or any portion of it. - - 11. Patents. - - A "contributor" is a copyright holder who authorizes use under this -License of the Program or a work on which the Program is based. The -work thus licensed is called the contributor's "contributor version". - - A contributor's "essential patent claims" are all patent claims -owned or controlled by the contributor, whether already acquired or -hereafter acquired, that would be infringed by some manner, permitted -by this License, of making, using, or selling its contributor version, -but do not include claims that would be infringed only as a -consequence of further modification of the contributor version. For -purposes of this definition, "control" includes the right to grant -patent sublicenses in a manner consistent with the requirements of -this License. - - Each contributor grants you a non-exclusive, worldwide, royalty-free -patent license under the contributor's essential patent claims, to -make, use, sell, offer for sale, import and otherwise run, modify and -propagate the contents of its contributor version. - - In the following three paragraphs, a "patent license" is any express -agreement or commitment, however denominated, not to enforce a patent -(such as an express permission to practice a patent or covenant not to -sue for patent infringement). To "grant" such a patent license to a -party means to make such an agreement or commitment not to enforce a -patent against the party. - - If you convey a covered work, knowingly relying on a patent license, -and the Corresponding Source of the work is not available for anyone -to copy, free of charge and under the terms of this License, through a -publicly available network server or other readily accessible means, -then you must either (1) cause the Corresponding Source to be so -available, or (2) arrange to deprive yourself of the benefit of the -patent license for this particular work, or (3) arrange, in a manner -consistent with the requirements of this License, to extend the patent -license to downstream recipients. "Knowingly relying" means you have -actual knowledge that, but for the patent license, your conveying the -covered work in a country, or your recipient's use of the covered work -in a country, would infringe one or more identifiable patents in that -country that you have reason to believe are valid. - - If, pursuant to or in connection with a single transaction or -arrangement, you convey, or propagate by procuring conveyance of, a -covered work, and grant a patent license to some of the parties -receiving the covered work authorizing them to use, propagate, modify -or convey a specific copy of the covered work, then the patent license -you grant is automatically extended to all recipients of the covered -work and works based on it. - - A patent license is "discriminatory" if it does not include within -the scope of its coverage, prohibits the exercise of, or is -conditioned on the non-exercise of one or more of the rights that are -specifically granted under this License. You may not convey a covered -work if you are a party to an arrangement with a third party that is -in the business of distributing software, under which you make payment -to the third party based on the extent of your activity of conveying -the work, and under which the third party grants, to any of the -parties who would receive the covered work from you, a discriminatory -patent license (a) in connection with copies of the covered work -conveyed by you (or copies made from those copies), or (b) primarily -for and in connection with specific products or compilations that -contain the covered work, unless you entered into that arrangement, -or that patent license was granted, prior to 28 March 2007. - - Nothing in this License shall be construed as excluding or limiting -any implied license or other defenses to infringement that may -otherwise be available to you under applicable patent law. - - 12. No Surrender of Others' Freedom. - - If conditions are imposed on you (whether by court order, agreement or -otherwise) that contradict the conditions of this License, they do not -excuse you from the conditions of this License. If you cannot convey a -covered work so as to satisfy simultaneously your obligations under this -License and any other pertinent obligations, then as a consequence you may -not convey it at all. For example, if you agree to terms that obligate you -to collect a royalty for further conveying from those to whom you convey -the Program, the only way you could satisfy both those terms and this -License would be to refrain entirely from conveying the Program. - - 13. Use with the GNU Affero General Public License. - - Notwithstanding any other provision of this License, you have -permission to link or combine any covered work with a work licensed -under version 3 of the GNU Affero General Public License into a single -combined work, and to convey the resulting work. The terms of this -License will continue to apply to the part which is the covered work, -but the special requirements of the GNU Affero General Public License, -section 13, concerning interaction through a network will apply to the -combination as such. - - 14. Revised Versions of this License. - - The Free Software Foundation may publish revised and/or new versions of -the GNU General Public License from time to time. Such new versions will -be similar in spirit to the present version, but may differ in detail to -address new problems or concerns. - - Each version is given a distinguishing version number. If the -Program specifies that a certain numbered version of the GNU General -Public License "or any later version" applies to it, you have the -option of following the terms and conditions either of that numbered -version or of any later version published by the Free Software -Foundation. If the Program does not specify a version number of the -GNU General Public License, you may choose any version ever published -by the Free Software Foundation. - - If the Program specifies that a proxy can decide which future -versions of the GNU General Public License can be used, that proxy's -public statement of acceptance of a version permanently authorizes you -to choose that version for the Program. - - Later license versions may give you additional or different -permissions. However, no additional obligations are imposed on any -author or copyright holder as a result of your choosing to follow a -later version. - - 15. Disclaimer of Warranty. - - THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY -APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT -HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY -OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, -THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR -PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM -IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF -ALL NECESSARY SERVICING, REPAIR OR CORRECTION. - - 16. Limitation of Liability. - - IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING -WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS -THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY -GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE -USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF -DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD -PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), -EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF -SUCH DAMAGES. - - 17. Interpretation of Sections 15 and 16. - - If the disclaimer of warranty and limitation of liability provided -above cannot be given local legal effect according to their terms, -reviewing courts shall apply local law that most closely approximates -an absolute waiver of all civil liability in connection with the -Program, unless a warranty or assumption of liability accompanies a -copy of the Program in return for a fee. - - END OF TERMS AND CONDITIONS - - How to Apply These Terms to Your New Programs - - If you develop a new program, and you want it to be of the greatest -possible use to the public, the best way to achieve this is to make it -free software which everyone can redistribute and change under these terms. - - To do so, attach the following notices to the program. It is safest -to attach them to the start of each source file to most effectively -state the exclusion of warranty; and each file should have at least -the "copyright" line and a pointer to where the full notice is found. - - fracture - Copyright (C) 2016 Jaron Kent-Dobias - - This program is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - This program is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with this program. If not, see <http://www.gnu.org/licenses/>. - -Also add information on how to contact you by electronic and paper mail. - - If the program does terminal interaction, make it output a short -notice like this when it starts in an interactive mode: - - fracture Copyright (C) 2016 Jaron Kent-Dobias - This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'. - This is free software, and you are welcome to redistribute it - under certain conditions; type `show c' for details. - -The hypothetical commands `show w' and `show c' should show the appropriate -parts of the General Public License. Of course, your program's commands -might be different; for a GUI interface, you would use an "about box". - - You should also get your employer (if you work as a programmer) or school, -if any, to sign a "copyright disclaimer" for the program, if necessary. -For more information on this, and how to apply and follow the GNU GPL, see -<http://www.gnu.org/licenses/>. - - The GNU General Public License does not permit incorporating your program -into proprietary programs. If your program is a subroutine library, you -may consider it more useful to permit linking proprietary applications with -the library. If this is what you want to do, use the GNU Lesser General -Public License instead of this License. But first, please read -<http://www.gnu.org/philosophy/why-not-lgpl.html>. diff --git a/README.md b/README.md deleted file mode 100644 index 2c1ec50..0000000 --- a/README.md +++ /dev/null @@ -1,4 +0,0 @@ -fracture -======== - -tools for simulating fracture in fuse networks
\ No newline at end of file diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt new file mode 100644 index 0000000..bf3c3b5 --- /dev/null +++ b/lib/CMakeLists.txt @@ -0,0 +1,4 @@ + +add_library(frac SHARED src/graph.cpp src/network.cpp) +target_include_directories(frac PUBLIC include) + diff --git a/lib/bound_set.c b/lib/bound_set.c deleted file mode 100644 index 32d0fb7..0000000 --- a/lib/bound_set.c +++ /dev/null @@ -1,110 +0,0 @@ - -#include "fracture.h" - -double th_p(double x, double y, double th) { - if (x >= 0 && y >= 0) - return th; - else if (x < 0 && y >= 0) - return M_PI - th; - else if (x < 0 && y < 0) - return th - M_PI; - else - return -th; -} - -double u_y(double x, double y) { - double r = sqrt(pow(x, 2) + pow(y, 2)); - double th = th_p(x, y, atan(fabs(y / x))); - - return sqrt(r) * sin(th / 2); -} - -void bound_set_embedded(double *bound, const graph_t *g, double notch_len) { - uint_t L = g->L; - - for (uint_t i = 0; i < L / 2; i++) { - double x1, y1, x2, y2, x3, y3, x4, y4; - x1 = (2. * i + 1.) / L - notch_len; - y1 = 0.5 - 1.; - x2 = (2. * i + 1.) / L - notch_len; - y2 = 0.5 - 0.; - y3 = (2. * i + 1.) / L - 0.5; - x3 = 0.5 - 1.; - y4 = (2. * i + 1.) / L - 0.5; - x4 = 0.5 - 0.; - - bound[g->b[g->bi[0] + i]] = u_y(x1, y1); - bound[g->b[g->bi[1] + i]] = u_y(x2, y2); - bound[g->b[g->bi[2] + i]] = u_y(x3, y3); - bound[g->b[g->bi[3] + i]] = u_y(x4, y4); - } -} - -bool is_in(uint_t len, uint_t *list, uint_t element) { - for (uint_t i = 0; i < len; i++) { - if (list[i] == element) { - return true; - } - } - return false; -} - -cholmod_dense *bound_set(const graph_t *g, bool vb, double notch_len, - cholmod_common *c) { - - uint_t dim = g->nv; - - if (vb && g->boundary != TORUS_BOUND) { - dim -= g->bi[g->nb]; - } else if (!vb) { - dim += 2; - } - - cholmod_dense *boundary = CHOL_F(zeros)(dim, 1, CHOLMOD_REAL, c); - double *bound = (double *)boundary->x; - - switch (g->boundary) { - case TORUS_BOUND: - for (uint_t i = 0; i < g->bi[1]; i++) { - uint_t be = g->b[i]; - uint_t v1 = g->ev[2 * be]; - uint_t v2 = g->ev[2 * be + 1]; - double v1y = g->vx[2 * v1 + 1]; - double v2y = g->vx[2 * v2 + 1]; - - uint_t ind = v1y < v2y ? 0 : 1; - - bound[g->ev[2 * be + ind]] += 1; - bound[g->ev[2 * be + !ind]] -= 1; - } - break; - /* -case EMBEDDED_BOUND: - bound_set_embedded(bound, g, notch_len); - break; - */ - default: - if (vb) { - for (uint_t i = 0; i < dim; i++) { - uint_t v = g->nbi[i]; - for (uint_t j = 0; j < g->vei[v + 1] - g->vei[v]; j++) { - uint_t e = g->ve[g->vei[v] + j]; - uint_t v0 = g->ev[2 * e]; - uint_t v1 = g->ev[2 * e + 1]; - - if (g->bq[v0] || g->bq[v1]) { - uint_t vv = v0 == v ? v1 : v0; - if (is_in(g->bi[1], g->b, vv)) { - bound[i]++; - } - } - } - } - } else { - bound[g->nv] = 1; - bound[g->nv + 1] = -1; - } - } - - return boundary; -} diff --git a/lib/break_edge.c b/lib/break_edge.c deleted file mode 100644 index e53f7c1..0000000 --- a/lib/break_edge.c +++ /dev/null @@ -1,50 +0,0 @@ - -#include "fracture.h" - -bool break_edge(net_t *net, uint_t edge, cholmod_common *c) { - net->fuses[edge] = true; - net->num_broken++; - - double *x = (double *)net->boundary_cond->x; - const graph_t *g = net->graph; - - uint_t v1 = net->graph->ev[2 * edge]; - uint_t v2 = net->graph->ev[2 * edge + 1]; - - uint_t s1 = v1 < v2 ? v1 : v2; - uint_t s2 = v1 < v2 ? v2 : v1; - - if (net->graph->boundary == TORUS_BOUND) { - if (net->graph->bq[edge]) { - double v1y = g->vx[2 * v1 + 1]; - double v2y = g->vx[2 * v2 + 1]; - uint_t ind = v1y < v2y ? 0 : 1; - x[g->ev[2 * edge + ind]] -= 1; - x[g->ev[2 * edge + !ind]] += 1; - } - if (net->factor != NULL) { - factor_update(net->factor, s1, s2, c); - } - } else if (net->voltage_bound) { - if (g->bq[v1] || g->bq[v2]) { - uint_t vv = g->bq[v1] ? v2 : v1; - uint_t uu = v1 == vv ? v2 : v1; - if (is_in(g->bi[1], g->b, uu)) { - x[g->bni[vv]] -= 1; - } - if (net->factor != NULL) { - factor_update2(net->factor, g->bni[vv], c); - } - } else { - if (net->factor != NULL) { - factor_update(net->factor, g->bni[s1], g->bni[s2], c); - } - } - } else { - if (net->factor != NULL) { - factor_update(net->factor, s1, s2, c); - } - } - - return true; -} diff --git a/lib/correlations.c b/lib/correlations.c deleted file mode 100644 index 98c5767..0000000 --- a/lib/correlations.c +++ /dev/null @@ -1,46 +0,0 @@ - -#include "fracture.h" - -double *get_corr(net_t *instance, uint_t **dists, cholmod_common *c) { - uint_t nv = instance->graph->dnv; - uint_t ne = instance->graph->ne; - uint_t *ev = instance->graph->dev; - bool nulldists = false; - if (dists == NULL) { - dists = graph_dists(instance->graph, false); - nulldists = true; - } - double *corr = calloc(nv, sizeof(double)); - uint_t *numat = calloc(nv, sizeof(uint_t)); - - for (uint_t i = 0; i < ne; i++) { - uint_t v1 = ev[2 * i]; - uint_t v2 = ev[2 * i + 1]; - for (uint_t j = 0; j < ne; j++) { - uint_t v3 = ev[2 * j]; - uint_t v4 = ev[2 * j + 1]; - uint_t dist1 = dists[v1][v3]; - uint_t dist2 = dists[v1][v4]; - uint_t dist3 = dists[v2][v3]; - uint_t dist4 = dists[v2][v4]; - uint_t dist = (dist1 + dist2 + dist3 + dist4) / 4; - corr[dist] += instance->fuses[i] && instance->fuses[j]; - numat[dist]++; - } - } - - for (uint_t i = 0; i < nv; i++) { - if (numat[i] > 0) { - corr[i] /= numat[i]; - } - } - - if (nulldists) { - for (int i = 0; i < nv; i++) { - free(dists[i]); - } - free(dists); - } - - return corr; -} diff --git a/lib/data.c b/lib/data.c deleted file mode 100644 index e39c450..0000000 --- a/lib/data.c +++ /dev/null @@ -1,35 +0,0 @@ - -#include "fracture.h" - -data_t *data_create(uint_t ne) { - data_t *data = malloc(1 * sizeof(data_t)); - assert(data != NULL); - - data->num_broken = 0; - - data->break_list = (uint_t *)malloc(ne * sizeof(uint_t)); - assert(data->break_list != NULL); - - data->extern_field = (long double *)malloc(ne * sizeof(long double)); - assert(data->extern_field != NULL); - - data->conductivity = (double *)malloc(ne * sizeof(double)); - assert(data->conductivity != NULL); - - return data; -} - -void data_free(data_t *data) { - free(data->break_list); - free(data->extern_field); - free(data->conductivity); - free(data); -} - -void data_update(data_t *data, uint_t last_broke, long double strength, - double conductivity) { - data->break_list[data->num_broken] = last_broke; - data->extern_field[data->num_broken] = strength; - data->conductivity[data->num_broken] = conductivity; - data->num_broken++; -} diff --git a/lib/factor_update.c b/lib/factor_update.c deleted file mode 100644 index ceaa01a..0000000 --- a/lib/factor_update.c +++ /dev/null @@ -1,69 +0,0 @@ - -#include "fracture.h" - -void factor_update(cholmod_factor *factor, uint_t v1, uint_t v2, - cholmod_common *c) { - uint_t n = factor->n; - - cholmod_sparse *update_mat = - CHOL_F(allocate_sparse)(n, n, 2, true, true, 0, CHOLMOD_REAL, c); - - uint_t s1, s2; - s1 = v1 < v2 ? v1 : v2; - s2 = v1 > v2 ? v1 : v2; - - int_t *pp = (int_t *)update_mat->p; - int_t *ii = (int_t *)update_mat->i; - double *xx = (double *)update_mat->x; - - for (uint_t i = 0; i <= s1; i++) { - pp[i] = 0; - } - - for (uint_t i = s1 + 1; i <= n; i++) { - pp[i] = 2; - } - - ii[0] = s1; - ii[1] = s2; - xx[0] = 1; - xx[1] = -1; - - cholmod_sparse *perm_update_mat = CHOL_F(submatrix)( - update_mat, factor->Perm, factor->n, NULL, -1, true, true, c); - - CHOL_F(updown)(false, perm_update_mat, factor, c); - - CHOL_F(free_sparse)(&perm_update_mat, c); - CHOL_F(free_sparse)(&update_mat, c); -} - -void factor_update2(cholmod_factor *factor, uint_t v, cholmod_common *c) { - uint_t n = factor->n; - - cholmod_sparse *update_mat = - CHOL_F(allocate_sparse)(n, n, 1, true, true, 0, CHOLMOD_REAL, c); - - int_t *pp = (int_t *)update_mat->p; - int_t *ii = (int_t *)update_mat->i; - double *xx = (double *)update_mat->x; - - for (uint_t i = 0; i <= v; i++) { - pp[i] = 0; - } - - for (uint_t i = v + 1; i <= n; i++) { - pp[i] = 1; - } - - ii[0] = v; - xx[0] = 1; - - cholmod_sparse *perm_update_mat = CHOL_F(submatrix)( - update_mat, factor->Perm, factor->n, NULL, -1, true, true, c); - - CHOL_F(updown)(false, perm_update_mat, factor, c); - - CHOL_F(free_sparse)(&perm_update_mat, c); - CHOL_F(free_sparse)(&update_mat, c); -} diff --git a/lib/fracture.h b/lib/fracture.h deleted file mode 100644 index 5eb0a1d..0000000 --- a/lib/fracture.h +++ /dev/null @@ -1,123 +0,0 @@ - -#pragma once - -#include <assert.h> -#include <cholmod.h> -#include <float.h> -#include <getopt.h> -#include <gsl/gsl_math.h> -#include <gsl/gsl_randist.h> -#include <gsl/gsl_rng.h> -#include <gsl/gsl_sf_erf.h> -#include <gsl/gsl_sf_exp.h> -#include <gsl/gsl_sf_log.h> -#include <inttypes.h> -#include <math.h> -#include <stdbool.h> -#include <stdint.h> -#include <stdio.h> -#include <stdlib.h> -#include <string.h> -#include <sys/types.h> -#include <unistd.h> - -#include <jst/graph.h> -#include <jst/rand.h> - -// these defs allow me to switch to long int cholmod in a sitch -#define int_t int -#define uint_t unsigned int -#define CINT_MAX INT_MAX -#define CHOL_F(x) cholmod_##x - -#define GSL_RAND_GEN gsl_rng_mt19937 - -typedef struct { - const graph_t *graph; - bool *fuses; - long double *thres; - double inf; - cholmod_dense *boundary_cond; - cholmod_factor *factor; - bool voltage_bound; - uint_t num_broken; - uint_t dim; - uint_t nep; - uint_t *evp; - cholmod_sparse *voltcurmat; -} net_t; - -typedef struct { - uint_t num_broken; - uint_t *break_list; - double *conductivity; - long double *extern_field; -} data_t; - -intptr_t *run_voronoi(uint_t num_coords, double *coords, bool periodic, - double xmin, double xmax, double ymin, double ymax); - -cholmod_sparse *gen_adjacency(const net_t *net, bool dual, bool use_gp, - bool symmetric, cholmod_common *c); - -cholmod_sparse *gen_laplacian(const net_t *net, cholmod_common *c); - -int edge_to_verts(uint_t width, bool periodic, uint_t edge, bool index); - -int dual_edge_to_verts(uint_t width, bool periodic, uint_t edge, bool index); - -double dual_vert_to_coord(uint_t width, bool periodic, uint_t vert, bool index); - -void factor_update(cholmod_factor *factor, uint_t v1, uint_t v2, - cholmod_common *c); -void factor_update2(cholmod_factor *factor, uint_t v1, cholmod_common *c); - -void net_notch(net_t *net, double notch_len, cholmod_common *c); -data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff); -double *net_voltages(const net_t *net, cholmod_common *c); -double *net_currents(const net_t *net, const double *voltages, - cholmod_common *c); -double net_conductivity(const net_t *net, const double *voltages); - -void update_boundary(net_t *instance, const double *avg_field); - -FILE *get_file(const char *prefix, uint_t width, uint_t crack, double beta, - uint_t iter, uint_t num_iter, uint_t num, bool read); - -double update_beta(double beta, uint_t width, const double *stress, - const double *damage, double bound_total); - -cholmod_sparse *gen_voltcurmat(uint_t num_edges, uint_t num_verts, - uint_t *edges_to_verts, cholmod_common *c); - -net_t *net_copy(const net_t *net, cholmod_common *c); - -void net_free(net_t *instance, cholmod_common *c); - -net_t *net_create(const graph_t *g, double inf, double beta, double notch_len, - bool vb, cholmod_common *c); - -bool break_edge(net_t *instance, uint_t edge, cholmod_common *c); - -components_t *get_clusters(net_t *instance); - -uint_t *get_cluster_dist(net_t *instance); - -void randfunc_flat(gsl_rng *r, double *x, double *y); -void randfunc_gaus(gsl_rng *r, double *x, double *y); - -double *get_corr(net_t *instance, uint_t **dists, cholmod_common *c); - -double *bin_values(graph_t *network, uint_t width, double *values); - -cholmod_dense *bound_set(const graph_t *g, bool vb, double notch_len, - cholmod_common *c); - -data_t *data_create(uint_t num_edges); -void data_free(data_t *data); -void data_update(data_t *data, uint_t last_broke, long double strength, - double conductivity); - -long double rand_dist_pow(const gsl_rng *r, double beta); - -bool is_in(uint_t len, uint_t *list, uint_t element); diff --git a/lib/gen_laplacian.c b/lib/gen_laplacian.c deleted file mode 100644 index 75dccce..0000000 --- a/lib/gen_laplacian.c +++ /dev/null @@ -1,142 +0,0 @@ - -#include "fracture.h" - -cholmod_sparse *gen_adjacency(const net_t *net, bool dual, bool use_gp, - bool symmetric, cholmod_common *c) { - const graph_t *g = net->graph; - uint_t nv; - uint_t ne; - uint_t nre; - uint_t *ev; - - if (use_gp) { - nv = net->dim; - ne = net->nep; - nre = (int_t)net->nep - (int_t)net->num_broken; - ev = net->evp; - } else if (dual) { - nv = g->dnv; - ne = g->ne; - nre = net->num_broken; - ev = g->dev; - } else { - nv = g->nv; - ne = g->ne; - nre = (int_t)g->ne - (int_t)net->num_broken; - ev = g->ev; - } - - uint_t nnz = nre; - - cholmod_triplet *t = - CHOL_F(allocate_triplet)(nv, nv, nnz, 1, CHOLMOD_REAL, c); - - int_t *ri = (int_t *)t->i; - int_t *ci = (int_t *)t->j; - double *ai = (double *)t->x; - - t->nnz = nnz; - - uint_t a = 0; - - for (uint_t i = 0; i < ne; i++) { - if ((net->fuses[i] && dual) || (!net->fuses[i] && !dual)) { - uint_t v1 = ev[2 * i]; - uint_t v2 = ev[2 * i + 1]; - - uint_t s1 = v1 < v2 ? v1 : v2; - uint_t s2 = v1 < v2 ? v2 : v1; - - ri[a] = s2; - ci[a] = s1; - ai[a] = 1; - - a++; - } - } - - cholmod_sparse *s = CHOL_F(triplet_to_sparse)(t, nnz, c); - CHOL_F(free_triplet)(&t, c); - - if (!symmetric) { - cholmod_sparse *tmp_s = CHOL_F(copy)(s, 0, 1, c); - CHOL_F(free_sparse)(&s, c); - s = tmp_s; - } - - return s; -} - -cholmod_sparse *gen_laplacian(const net_t *net, cholmod_common *c) { - const graph_t *g = net->graph; - uint_t nv = net->dim; - uint_t ne = net->nep; - uint_t *ev = net->evp; - - uint_t nnz = nv; - - cholmod_triplet *temp_m = - CHOL_F(allocate_triplet)(nv, nv, nnz, 1, CHOLMOD_REAL, c); - - int_t *rowind = (int_t *)temp_m->i; - int_t *colind = (int_t *)temp_m->j; - double *acoo = (double *)temp_m->x; - - temp_m->nnz = nnz; - - for (uint_t i = 0; i < nv; i++) { - rowind[i] = i; - colind[i] = i; - acoo[i] = 0; - } - - cholmod_sparse *adjacency = gen_adjacency(net, false, true, true, c); - - for (uint_t i = 0; i < ne; i++) { - if (!net->fuses[i]) { - uint_t v1 = ev[2 * i]; - uint_t v2 = ev[2 * i + 1]; - - acoo[v1]++; - acoo[v2]++; - } - } - - if (net->voltage_bound && g->boundary != TORUS_BOUND) { - for (uint_t i = 0; i < net->dim; i++) { - uint_t v = g->nbi[i]; - for (uint_t j = 0; j < g->vei[v + 1] - g->vei[v]; j++) { - uint_t e = g->ve[g->vei[v] + j]; - uint_t v0 = g->ev[2 * e]; - uint_t v1 = g->ev[2 * e + 1]; - - if (g->bq[v0] || g->bq[v1]) { - acoo[i]++; - } - } - } - } else { - acoo[0]++; - } - - for (uint_t i = 0; i < nv; i++) { - if (acoo[i] == 0) - acoo[i]++; - } - - // assert(CHOL_F(check_triplet)(temp_m, c)); - - cholmod_sparse *t_out = CHOL_F(triplet_to_sparse)(temp_m, nnz, c); - // assert(CHOL_F(check_sparse)(t_out, c)); - - double alpha[2] = {1, 0}; - double beta[2] = {-1, 0}; - cholmod_sparse *laplacian = - CHOL_F(add)(t_out, adjacency, alpha, beta, 1, 1, c); - - CHOL_F(free_sparse)(&t_out, c); - CHOL_F(free_sparse)(&adjacency, c); - CHOL_F(free_triplet)(&temp_m, c); - - return laplacian; -} diff --git a/lib/gen_voltcurmat.c b/lib/gen_voltcurmat.c deleted file mode 100644 index fede836..0000000 --- a/lib/gen_voltcurmat.c +++ /dev/null @@ -1,36 +0,0 @@ - -#include "fracture.h" - -cholmod_sparse *gen_voltcurmat(unsigned int num_edges, unsigned int num_verts, - unsigned int *edges_to_verts, - cholmod_common *c) { - - cholmod_triplet *t_m = CHOL_F(allocate_triplet)( - num_edges, num_verts, num_edges * 2, 0, CHOLMOD_REAL, c); - assert(t_m != NULL); - - int_t *rowind = (int_t *)t_m->i; - int_t *colind = (int_t *)t_m->j; - double *acoo = (double *)t_m->x; - - for (unsigned int i = 0; i < num_edges; i++) { - unsigned int v1 = edges_to_verts[2 * i]; - unsigned int v2 = edges_to_verts[2 * i + 1]; - rowind[2 * i] = i; - rowind[2 * i + 1] = i; - colind[2 * i] = v1; - colind[2 * i + 1] = v2; - acoo[2 * i] = 1; - acoo[2 * i + 1] = -1; - } - - t_m->nnz = num_edges * 2; - - assert(CHOL_F(check_triplet)(t_m, c)); - - cholmod_sparse *m = CHOL_F(triplet_to_sparse)(t_m, num_edges * 2, c); - - CHOL_F(free_triplet)(&t_m, c); - - return m; -} diff --git a/lib/geometry.c b/lib/geometry.c deleted file mode 100644 index 2c1987b..0000000 --- a/lib/geometry.c +++ /dev/null @@ -1,55 +0,0 @@ - -#include "fracture.h" - -int edge_to_verts(unsigned int width, bool periodic, unsigned int edge, - bool index) { - assert(edge < pow(width, 2)); - - int x = edge / width + 1; - int y = edge % width + 1; - - if (periodic) { - return (((index ^ (x % 2)) + 2 * ((y + (index ^ (!(x % 2)))) / 2) - 1) % - width + - (x - index) * width) / - 2; - } else { - return ((index ^ (x % 2)) + 2 * ((y + (index ^ (!(x % 2)))) / 2) + - (x - index) * (width + 1) - 1) / - 2; - } -} - -int dual_edge_to_verts(unsigned int width, bool periodic, unsigned int edge, - bool index) { - assert(edge < pow(width, 2)); - - int x = edge / width + 1; - int y = edge % width + 1; - - if (periodic) { - return (((index ^ (!(x % 2))) + 2 * ((y + (index ^ (x % 2))) / 2) - 1) % - width + - (x - index) * width) / - 2; - } else { - return ((index ^ (!(x % 2))) + 2 * ((y + (index ^ (x % 2))) / 2) + - (x - index) * (width + 1) - 1) / - 2; - } -} - -double dual_vert_to_coord(unsigned int width, bool periodic, unsigned int vert, - bool index) { - if (periodic) { - if (index) - return (2 * vert) % width + (2 * vert / width) % 2; - else - return 2 * vert / width; - } else { - if (index) - return (2 * vert) % (width + 1); - else - return (2 * vert) / (width + 1); - } -} diff --git a/lib/get_dual_clusters.c b/lib/get_dual_clusters.c deleted file mode 100644 index 9336106..0000000 --- a/lib/get_dual_clusters.c +++ /dev/null @@ -1,31 +0,0 @@ - -#include "fracture.h" - -components_t *get_clusters(net_t *instance) { - components_t *c = - graph_components_get(instance->graph, instance->fuses, true); - return c; -} - -unsigned int *get_cluster_dist(net_t *instance) { - components_t *c = get_clusters(instance); - unsigned int *cluster_dist = - (unsigned int *)calloc(instance->graph->dnv, sizeof(unsigned int)); - - for (uint32_t i = 1; i <= c->n; i++) { - unsigned int num_in_cluster = 0; - for (unsigned int j = 0; j < instance->graph->dnv; j++) { - if (c->labels[j] == i) - num_in_cluster++; - } - - if (num_in_cluster == 0) - break; - - cluster_dist[num_in_cluster - 1]++; - } - - graph_components_free(c); - - return cluster_dist; -} diff --git a/lib/include/graph.hpp b/lib/include/graph.hpp new file mode 100644 index 0000000..80cdd66 --- /dev/null +++ b/lib/include/graph.hpp @@ -0,0 +1,29 @@ + +#pragma once + +#include <cmath> +#include <vector> +#include <array> + +class graph { + public: + typedef struct coordinate { + double x; + double y; + } coordinate; + + typedef struct vertex { + coordinate r; + } vertex; + + typedef std::array<unsigned int, 2> edge; + + std::vector<vertex> vertices; + std::vector<edge> edges; + + std::vector<vertex> dual_vertices; + std::vector<edge> dual_edges; + + graph(unsigned int L); +}; + diff --git a/lib/include/hooks.hpp b/lib/include/hooks.hpp new file mode 100644 index 0000000..350f318 --- /dev/null +++ b/lib/include/hooks.hpp @@ -0,0 +1,12 @@ + +#pragma once + +class network; + +class hooks { + public: + virtual void pre_fracture(const network&) {}; + virtual void bond_broken(const network&, const std::pair<double, std::vector<double>>&, unsigned int) {}; + virtual void post_fracture(network&) {}; // post fracture hook can be destructive +}; + diff --git a/lib/include/network.hpp b/lib/include/network.hpp new file mode 100644 index 0000000..abf88cd --- /dev/null +++ b/lib/include/network.hpp @@ -0,0 +1,48 @@ + +#pragma once + +#include <vector> +#include <functional> +#include <utility> +#include <random> + +#include "cholmod.h" + +#include "graph.hpp" +#include "hooks.hpp" + +#ifdef FRACTURE_LONGINT + +#define CHOL_F(x) cholmod_l_##x +#define CHOL_INT long int + +#else + +#define CHOL_F(x) cholmod_##x +#define CHOL_INT int + +#endif + +class network { + private: + cholmod_dense *b; + cholmod_factor *factor; + cholmod_sparse *voltcurmat; + cholmod_common *c; + + public: + const graph& G; + std::vector<bool> fuses; + std::vector<long double> thresholds; + + network(const graph&, cholmod_common*); + network(const network &other); + ~network(); + + void set_thresholds(double, std::mt19937&); + void break_edge(unsigned int); + void add_edge(unsigned int); + std::pair<double, std::vector<double>> currents(); + void fracture(hooks&, double cutoff = 1e-13); +}; + diff --git a/lib/net.c b/lib/net.c deleted file mode 100644 index ac89415..0000000 --- a/lib/net.c +++ /dev/null @@ -1,146 +0,0 @@ - -#include "fracture.h" - -long double *get_thres(uint_t ne, double beta) { - long double *thres = (long double *)malloc(ne * sizeof(long double)); - assert(thres != NULL); - - gsl_rng *r = gsl_rng_alloc(GSL_RAND_GEN); - gsl_rng_set(r, jst_rand_seed()); - - for (uint_t i = 0; i < ne; i++) { - thres[i] = rand_dist_pow(r, beta); - } - - gsl_rng_free(r); - - return thres; -} - -void net_notch(net_t *net, double notch_len, cholmod_common *c) { - for (uint_t i = 0; i < net->graph->ne; i++) { - uint_t v1, v2; - double v1x, v1y, v2x, v2y, dy; - bool crosses_center, not_wrapping, correct_length; - - v1 = net->graph->ev[2 * i]; - v2 = net->graph->ev[2 * i + 1]; - - v1x = net->graph->vx[2 * v1]; - v1y = net->graph->vx[2 * v1 + 1]; - v2x = net->graph->vx[2 * v2]; - v2y = net->graph->vx[2 * v2 + 1]; - - dy = v1y - v2y; - - crosses_center = (v1y >= 0.5 && v2y <= 0.5) || (v1y <= 0.5 && v2y >= 0.5); - not_wrapping = fabs(dy) < 0.5; - // correct_length = v1x + dx / dy * (v1y - 0.5) <= notch_len; - correct_length = v1x < notch_len && v2x < notch_len; - - if (crosses_center && not_wrapping && correct_length) { - break_edge(net, i, c); - } - } -} - -net_t *net_create(const graph_t *g, double inf, double beta, double notch_len, - bool vb, cholmod_common *c) { - net_t *net = (net_t *)calloc(1, sizeof(net_t)); - assert(net != NULL); - - net->graph = g; - net->num_broken = 0; - - net->fuses = (bool *)calloc(g->ne, sizeof(bool)); - assert(net->fuses != NULL); - - net->thres = get_thres(g->ne, beta); - net->inf = inf; - - net->dim = g->nv; - - if (g->boundary == TORUS_BOUND) { - net->nep = g->ne; - net->evp = (uint_t *)malloc(2 * g->ne * sizeof(uint_t)); - memcpy(net->evp, g->ev, 2 * g->ne * sizeof(uint_t)); - } else { - if (vb) { - net->dim -= g->bi[g->nb]; - net->evp = (uint_t *)malloc(2 * g->ne * sizeof(uint_t)); - net->nep = 0; - for (uint_t i = 0; i < g->ne; i++) { - if (!(g->bq[g->ev[2 * i]] || g->bq[g->ev[2 * i + 1]])) { - net->evp[2 * net->nep] = g->bni[g->ev[2 * i]]; - net->evp[2 * net->nep + 1] = g->bni[g->ev[2 * i + 1]]; - net->nep++; - } - } - } else { - net->dim += 2; - net->evp = (uint_t *)malloc(2 * (g->ne + g->bi[2]) * sizeof(uint_t)); - memcpy(net->evp, g->ev, 2 * g->ne * sizeof(uint_t)); - net->nep = g->ne + g->bi[2]; - - for (uint_t i = 0; i < 2; i++) { - for (uint_t j = 0; j < g->bi[i + 1] - g->bi[i]; j++) { - net->evp[2 * (g->ne + g->bi[i] + j)] = g->b[g->bi[i] + j]; - net->evp[2 * (g->ne + g->bi[i] + j) + 1] = g->nv + i; - } - } - } - } - - net->voltage_bound = vb; - net->boundary_cond = bound_set(g, vb, notch_len, c); - - net_notch(net, notch_len, c); - - { - cholmod_sparse *laplacian = gen_laplacian(net, c); - net->factor = CHOL_F(analyze)(laplacian, c); - CHOL_F(factorize)(laplacian, net->factor, c); - CHOL_F(free_sparse)(&laplacian, c); - } - - net->voltcurmat = gen_voltcurmat(g->ne, g->nv, g->ev, c); - - return net; -} - -net_t *net_copy(const net_t *net, cholmod_common *c) { - net_t *net_copy = (net_t *)calloc(1, sizeof(net_t)); - assert(net_copy != NULL); - memcpy(net_copy, net, sizeof(net_t)); - - size_t fuses_size = (net->graph)->ne * sizeof(bool); - net_copy->fuses = (bool *)malloc(fuses_size); - assert(net_copy->fuses != NULL); - memcpy(net_copy->fuses, net->fuses, fuses_size); - - size_t thres_size = (net->graph)->ne * sizeof(long double); - net_copy->thres = (long double *)malloc(thres_size); - assert(net_copy->thres != NULL); - memcpy(net_copy->thres, net->thres, thres_size); - - size_t evp_size = 2 * net->nep * sizeof(uint_t); - net_copy->evp = (uint_t *)malloc(thres_size); - assert(net_copy->evp != NULL); - memcpy(net_copy->evp, net->evp, evp_size); - - net_copy->boundary_cond = CHOL_F(copy_dense)(net->boundary_cond, c); - net_copy->factor = CHOL_F(copy_factor)(net->factor, c); - net_copy->voltcurmat = CHOL_F(copy_sparse)(net->voltcurmat, c); - - return net_copy; -} - -void net_free(net_t *net, cholmod_common *c) { - free(net->fuses); - free(net->thres); - CHOL_F(free_dense)(&(net->boundary_cond), c); - CHOL_F(free_factor)(&(net->factor), c); - CHOL_F(free_sparse)(&(net->voltcurmat), c); - free(net->evp); - free(net); -} diff --git a/lib/net_conductivity.c b/lib/net_conductivity.c deleted file mode 100644 index 61148da..0000000 --- a/lib/net_conductivity.c +++ /dev/null @@ -1,35 +0,0 @@ - -#include "fracture.h" - -double net_conductivity(const net_t *net, const double *voltages) { - if (net->voltage_bound) { - // the voltage drop across the network is fixed to one with voltage - // boundary conditions, so the conductivity is the total current flowing - double tot_cur = 0; - for (uint_t i = 0; i < net->graph->num_spanning_edges; i++) { - uint_t e = net->graph->spanning_edges[i]; - - if (!net->fuses[e]) { - uint_t v1, v2, s1, s2; - double v1y, v2y; - - v1 = net->graph->ev[2 * e]; - v2 = net->graph->ev[2 * e + 1]; - - v1y = net->graph->vx[2 * v1 + 1]; - v2y = net->graph->vx[2 * v2 + 1]; - - s1 = v1y < v2y ? v1 : v2; - s2 = v1y < v2y ? v2 : v1; - - tot_cur += voltages[s1] - voltages[s2]; - } - } - - return fabs(tot_cur); - } else { - // the current across the network is fixed to one with current boundary - // conditions, so the conductivity is the inverse of the total voltage drop - return 1 / fabs(voltages[net->graph->nv] - voltages[net->graph->nv + 1]); - } -} diff --git a/lib/net_currents.c b/lib/net_currents.c deleted file mode 100644 index 32dba04..0000000 --- a/lib/net_currents.c +++ /dev/null @@ -1,52 +0,0 @@ - -#include "fracture.h" - -double *net_currents(const net_t *net, const double *voltages, - cholmod_common *c) { - uint_t ne = net->graph->ne; - uint_t dim = net->graph->nv; - cholmod_sparse *voltcurmat = net->voltcurmat; - - cholmod_dense *x = CHOL_F(allocate_dense)(dim, 1, dim, CHOLMOD_REAL, c); - - double *tmp_x = x->x; - x->x = (void *)voltages; - - cholmod_dense *y = CHOL_F(allocate_dense)(ne, 1, ne, CHOLMOD_REAL, c); - - double alpha[2] = {1, 0}; - double beta[2] = {0, 0}; - CHOL_F(sdmult)(voltcurmat, 0, alpha, beta, x, y, c); - - double *currents = (double *)malloc(ne * sizeof(double)); - - for (int i = 0; i < ne; i++) { - if (net->fuses[i]) { - currents[i] = 0; - } else { - currents[i] = ((double *)y->x)[i]; - } - } - - if (net->graph->boundary == TORUS_BOUND) { - for (uint_t i = 0; i < net->graph->bi[1]; i++) { - uint_t e = net->graph->b[i]; - uint_t v1 = net->graph->ev[2 * e]; - uint_t v2 = net->graph->ev[2 * e + 1]; - double v1y = net->graph->vx[2 * v1 + 1]; - double v2y = net->graph->vx[2 * v2 + 1]; - - if (v1y > v2y) { - currents[e] += 1; - } else { - currents[e] -= 1; - } - } - } - - x->x = tmp_x; - CHOL_F(free_dense)(&x, c); - CHOL_F(free_dense)(&y, c); - - return currents; -} diff --git a/lib/net_fracture.c b/lib/net_fracture.c deleted file mode 100644 index 65ede9b..0000000 --- a/lib/net_fracture.c +++ /dev/null @@ -1,68 +0,0 @@ - -#include "fracture.h" - -uint_t get_next_broken(net_t *net, double *currents, double cutoff) { - uint_t max_pos = UINT_MAX; - long double max_val = 0; - - for (uint_t i = 0; i < net->graph->ne; i++) { - long double current = fabs(currents[i]); - bool broken = net->fuses[i]; - - if (!broken && current > cutoff) { - long double val = current / net->thres[i]; - - if (val > max_val) { - max_val = val; - max_pos = i; - } - } - } - - if (max_pos == UINT_MAX) { - printf("GET_NEXT_BROKEN: currents is zero or NaN, no max_val found\n"); - exit(EXIT_FAILURE); - } - - return max_pos; -} - -data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff) { - data_t *data = data_create(net->graph->ne); - - uint_t n = 0; - while (true) { - n++; - double *voltages = net_voltages(net, c); - double *currents = net_currents(net, voltages, c); - - double conductivity = net_conductivity(net, voltages); - - if (conductivity < cutoff) { - free(voltages); - free(currents); - break; - } - - uint_t last_broke = get_next_broken(net, currents, cutoff); - - long double sim_current; - - if (net->voltage_bound) { - sim_current = conductivity; - } else { - sim_current = 1; - } - - data_update(data, last_broke, fabsl(sim_current * (net->thres)[last_broke] / - ((long double)currents[last_broke])), - conductivity); - - free(voltages); - free(currents); - - break_edge(net, last_broke, c); - } - - return data; -} diff --git a/lib/net_voltages.c b/lib/net_voltages.c deleted file mode 100644 index d456a65..0000000 --- a/lib/net_voltages.c +++ /dev/null @@ -1,39 +0,0 @@ - -#include "fracture.h" - -double *net_voltages(const net_t *net, cholmod_common *c) { - cholmod_dense *b = net->boundary_cond; - cholmod_factor *factor = net->factor; - - cholmod_dense *x = CHOL_F(solve)(CHOLMOD_A, factor, b, c); - - if (((double *)x->x)[0] != ((double *)x->x)[0]) { - printf("GET_VOLTAGE: value is NaN\n"); - exit(EXIT_FAILURE); - } - - double *t_voltages = (double *)x->x; - x->x = NULL; - CHOL_F(free_dense)(&x, c); - - const graph_t *g = net->graph; - - if (g->boundary == TORUS_BOUND) { - return t_voltages; - } else if (net->voltage_bound) { - double *voltages = (double *)malloc(g->nv * sizeof(double)); - for (uint_t i = 0; i < g->nv - g->bi[g->nb]; i++) { - voltages[g->nbi[i]] = t_voltages[i]; - } - for (uint_t i = 0; i < 2; i++) { - for (uint_t j = 0; j < g->bi[i + 1] - g->bi[i]; j++) { - voltages[g->b[g->bi[i] + j]] = 1 - i; - } - } - - free(t_voltages); - return voltages; - } else { - return t_voltages; - } -} diff --git a/lib/rand.c b/lib/rand.c deleted file mode 100644 index 8ba2b2b..0000000 --- a/lib/rand.c +++ /dev/null @@ -1,15 +0,0 @@ - -#include "fracture.h" - -long double rand_dist_pow(const gsl_rng *r, double beta) { - long double x = 0; - - // underflow means that for very small beta x is sometimes identically zero, - // which causes problems - while (x == 0.0) { - long double y = logl(gsl_rng_uniform_pos(r)) / beta; - x = expl(y); - } - - return x; -} diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp new file mode 100644 index 0000000..e76947b --- /dev/null +++ b/lib/src/graph.cpp @@ -0,0 +1,38 @@ + +#include "graph.hpp" + +graph::graph(unsigned int L) { + double dx = 1.0 / L; + unsigned int nv = 2 * pow(L / 2, 2); + unsigned int ne = pow(L, 2); + + vertices.resize(nv); + edges.resize(ne); + + dual_vertices.resize(nv); + dual_edges.resize(ne); + + for (unsigned int i = 0; i < nv; i++) { + vertices[i].r.x = dx * ((1 + i / (L / 2)) % 2 + 2 * (i % (L / 2))); + vertices[i].r.y = dx * (i / (L / 2)); + + dual_vertices[i].r.x = dx * ((i / (L / 2)) % 2 + 2 * (i % (L / 2))); + dual_vertices[i].r.y = dx * (i / (L / 2)); + } + + for (unsigned int i = 0; i < ne; i++) { + unsigned int x = i / L; + unsigned int y = i % L; + + unsigned int v1 = (L * x) / 2 + ((y + x % 2) / 2) % (L / 2); + unsigned int v2 = ((L * (x + 1)) / 2 + ((y + (x + 1) % 2) / 2) % (L / 2)) % nv; + + edges[i] = {v1, v2}; + + unsigned int dv1 = (L * x) / 2 + ((y + (x + 1) % 2) / 2) % (L / 2); + unsigned int dv2 = ((L * (x + 1)) / 2 + ((y + x % 2) / 2) % (L / 2)) % nv; + + dual_edges[i] = {dv1, dv2}; + } +} + diff --git a/lib/src/network.cpp b/lib/src/network.cpp new file mode 100644 index 0000000..8af140b --- /dev/null +++ b/lib/src/network.cpp @@ -0,0 +1,280 @@ + +#include "network.hpp" + +network::network(const graph& G, cholmod_common *c) : G(G), c(c), fuses(G.edges.size(), false), thresholds(G.edges.size(), 1) { + b = CHOL_F(zeros)(G.vertices.size(), 1, CHOLMOD_REAL, c); + for (unsigned int i = 0; i < G.edges.size(); i++) { + double v0y = G.vertices[G.edges[i][0]].r.y; + double v1y = G.vertices[G.edges[i][1]].r.y; + + if (fabs(v0y - v1y) > 0.5) { + bool ind = v0y < v1y ? 0 : 1; + + ((double *)b->x)[G.edges[i][ind]] += 1.0; + ((double *)b->x)[G.edges[i][!ind]] -= 1.0; + } + } + + unsigned int nnz = G.vertices.size() + G.edges.size(); + + cholmod_triplet *t = CHOL_F(allocate_triplet)(G.vertices.size(), G.vertices.size(), nnz, 1, CHOLMOD_REAL, c); + + t->nnz = nnz; + + for (unsigned int i = 0; i < G.vertices.size(); i++) { + ((CHOL_INT *)t->i)[i] = i; + ((CHOL_INT *)t->j)[i] = i; + ((double *)t->x)[i] = 0.0; + } + + unsigned int terms = G.vertices.size(); + + for (unsigned int i = 0; i < G.edges.size(); i++) { + unsigned int v0 = G.edges[i][0]; + unsigned int v1 = G.edges[i][1]; + + unsigned int s0 = v0 < v1 ? v0 : v1; + unsigned int s1 = v0 < v1 ? v1 : v0; + + ((CHOL_INT *)t->i)[terms] = v1; + ((CHOL_INT *)t->j)[terms] = v0; + ((double *)t->x)[terms] = -1.0; + + ((double *)t->x)[v0]++; + ((double *)t->x)[v1]++; + + terms++; + } + + ((double *)t->x)[0]++; + + cholmod_sparse *laplacian = CHOL_F(triplet_to_sparse)(t, nnz, c); + CHOL_F(free_triplet)(&t, c); + factor = CHOL_F(analyze)(laplacian, c); + CHOL_F(factorize)(laplacian, factor, c); + CHOL_F(free_sparse)(&laplacian, c); + + t = CHOL_F(allocate_triplet)(G.edges.size(), G.vertices.size(), 2 * G.edges.size(), 0, CHOLMOD_REAL, c); + + + t->nnz = 2 * G.edges.size(); + + for (unsigned int i = 0; i < G.edges.size(); i++) { + ((CHOL_INT *)t->i)[2 * i] = i; + ((CHOL_INT *)t->j)[2 * i] = G.edges[i][0]; + ((double *)t->x)[2 * i] = 1.0; + + ((CHOL_INT *)t->i)[2 * i + 1] = i; + ((CHOL_INT *)t->j)[2 * i + 1] = G.edges[i][1]; + ((double *)t->x)[2 * i + 1] = -1.0; + } + + voltcurmat = CHOL_F(triplet_to_sparse)(t, 2 * G.edges.size(), c); + + CHOL_F(free_triplet)(&t, c); +} + +network::network(const network &other) : G(other.G), thresholds(other.thresholds), fuses(other.fuses), c(other.c) { + b = CHOL_F(copy_dense)(other.b, c); + factor = CHOL_F(copy_factor)(other.factor, c); + voltcurmat = CHOL_F(copy_sparse)(other.voltcurmat, c); +} + +network::~network() { + CHOL_F(free_sparse)(&voltcurmat, c); + CHOL_F(free_dense)(&b, c); + CHOL_F(free_factor)(&factor, c); +} + +void network::set_thresholds(double beta, std::mt19937& rng) { + std::uniform_real_distribution<long double> d(0.0, 1.0); + + for (unsigned int i = 0; i < G.edges.size(); i++) { + long double x = 0.0; + + while (x == 0.0) { + x = expl(logl(d(rng)) / (long double)beta); + } + + thresholds[i] = x; + } +} + +void network::add_edge(unsigned int e) { + fuses[e] = false; + unsigned int v0 = G.edges[e][0]; + unsigned int v1 = G.edges[e][1]; + + unsigned int n = factor->n; + + cholmod_sparse *update_mat = CHOL_F(allocate_sparse)(n, n, 2, true, true, 0, CHOLMOD_REAL, c); + + unsigned int s1, s2; + s1 = v0 < v1 ? v0 : v1; + s2 = v0 < v1 ? v1 : v0; + + CHOL_INT *pp = (CHOL_INT *)update_mat->p; + CHOL_INT *ii = (CHOL_INT *)update_mat->i; + double *xx = (double *)update_mat->x; + + for (unsigned int i = 0; i <= s1; i++) { + pp[i] = 0; + } + + for (unsigned int i = s1 + 1; i <= n; i++) { + pp[i] = 2; + } + + ii[0] = s1; + ii[1] = s2; + xx[0] = -1; + xx[1] = 1; + + cholmod_sparse *perm_update_mat = CHOL_F(submatrix)( + update_mat, (CHOL_INT *)factor->Perm, factor->n, NULL, -1, true, true, c); + + CHOL_F(updown)(false, perm_update_mat, factor, c); + + CHOL_F(free_sparse)(&perm_update_mat, c); + CHOL_F(free_sparse)(&update_mat, c); + + double v0y = G.vertices[v0].r.y; + double v1y = G.vertices[v1].r.y; + + if (fabs(v0y - v1y) > 0.5) { + bool ind = v0y < v1y ? 0 : 1; + + ((double *)b->x)[G.edges[e][ind]] += 1.0; + ((double *)b->x)[G.edges[e][!ind]] -= 1.0; + } +} + +void network::break_edge(unsigned int e) { + fuses[e] = true; + unsigned int v0 = G.edges[e][0]; + unsigned int v1 = G.edges[e][1]; + + unsigned int n = factor->n; + + cholmod_sparse *update_mat = CHOL_F(allocate_sparse)(n, n, 2, true, true, 0, CHOLMOD_REAL, c); + + unsigned int s1, s2; + s1 = v0 < v1 ? v0 : v1; + s2 = v0 < v1 ? v1 : v0; + + CHOL_INT *pp = (CHOL_INT *)update_mat->p; + CHOL_INT *ii = (CHOL_INT *)update_mat->i; + double *xx = (double *)update_mat->x; + + for (unsigned int i = 0; i <= s1; i++) { + pp[i] = 0; + } + + for (unsigned int i = s1 + 1; i <= n; i++) { + pp[i] = 2; + } + + ii[0] = s1; + ii[1] = s2; + xx[0] = 1; + xx[1] = -1; + + cholmod_sparse *perm_update_mat = CHOL_F(submatrix)( + update_mat, (CHOL_INT *)factor->Perm, factor->n, NULL, -1, true, true, c); + + CHOL_F(updown)(false, perm_update_mat, factor, c); + + CHOL_F(free_sparse)(&perm_update_mat, c); + CHOL_F(free_sparse)(&update_mat, c); + + double v0y = G.vertices[v0].r.y; + double v1y = G.vertices[v1].r.y; + + if (fabs(v0y - v1y) > 0.5) { + bool ind = v0y < v1y ? 0 : 1; + + ((double *)b->x)[G.edges[e][ind]] -= 1.0; + ((double *)b->x)[G.edges[e][!ind]] += 1.0; + } +} + +std::pair<double, std::vector<double>> network::currents() { + cholmod_dense *x = CHOL_F(solve)(CHOLMOD_A, factor, b, c); + + if (((double *)x->x)[0] != ((double *)x->x)[0]) { + exit(2); + } + + cholmod_dense *y = CHOL_F(allocate_dense)(G.edges.size(), 1, G.edges.size(), CHOLMOD_REAL, c); + + double alpha[2] = {1, 0}; + double beta[2] = {0, 0}; + CHOL_F(sdmult)(voltcurmat, 0, alpha, beta, x, y, c); + + std::vector<double> currents(G.edges.size()); + + double total_current = 0; + + for (int i = 0; i < G.edges.size(); i++) { + if (fuses[i]) { + currents[i] = 0; + } else { + currents[i] = ((double *)y->x)[i]; + + double v0y = G.vertices[G.edges[i][0]].r.y; + double v1y = G.vertices[G.edges[i][1]].r.y; + + if (fabs(v0y - v1y) > 0.5) { + if (v0y > v1y) { + currents[i] += 1.0; + total_current += currents[i]; + } else { + currents[i] -= 1.0; + total_current -= currents[i]; + } + } + } + + } + + CHOL_F(free_dense)(&x, c); + CHOL_F(free_dense)(&y, c); + + return std::make_pair(total_current, currents); +} + +void network::fracture(hooks& m, double cutoff) { + m.pre_fracture(*this); + + while (true) { + auto currents = this->currents(); + + if (currents.first < cutoff) { + break; + } + + unsigned int max_pos = UINT_MAX; + long double max_val = 0; + + for (unsigned int i = 0; i < G.edges.size(); i++) { + if (!fuses[i] && fabs(currents.second[i]) > cutoff) { + long double val = (long double)fabs(currents.second[i]) / thresholds[i]; + if (val > max_val) { + max_val = val; + max_pos = i; + } + } + } + + if (max_pos == UINT_MAX) { + exit(3); + } + + this->break_edge(max_pos); + + m.bond_broken(*this, currents, max_pos); + } + + m.post_fracture(*this); +} + diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt new file mode 100644 index 0000000..57a817b --- /dev/null +++ b/src/CMakeLists.txt @@ -0,0 +1,7 @@ + +add_library(measurements measurements.cpp) +target_link_libraries(measurements frac) + +add_executable(fracture fracture.cpp) +target_link_libraries(fracture frac measurements fftw3 cholmod) + diff --git a/src/anal_process.c b/src/anal_process.c deleted file mode 100644 index de27571..0000000 --- a/src/anal_process.c +++ /dev/null @@ -1,135 +0,0 @@ - -#include "fracture.h" -#include <gsl/gsl_blas.h> -#include <gsl/gsl_matrix.h> -#include <gsl/gsl_sf_erf.h> -#include <gsl/gsl_sf_laguerre.h> -#include <gsl/gsl_vector.h> - -void mom(uint_t len, uint_t n, uint32_t *data, double result[2]) { - uint_t total = 0; - double mom = 0; - double mom_err = 0; - - for (uint_t i = 0; i < len; i++) { - uint32_t datai = data[i]; - double in = pow(i, n); - - total += datai; - mom += in * datai; - mom_err += gsl_pow_2(in) * datai; - } - - double momf = mom / total; - double momf_err = momf * sqrt(mom_err / gsl_pow_2(mom) + 1 / total); - - result[0] = momf; - result[1] = momf_err; -} - -int main(int argc, char *argv[]) { - uint_t nc = argc - 1; - uint_t *Ls_c = (uint_t *)malloc(nc * sizeof(uint_t)); - double *betas_c = (double *)malloc(nc * sizeof(double)); - double *vals_c1 = (double *)malloc(nc * sizeof(double)); - double *errors_c1 = (double *)malloc(nc * sizeof(double)); - double *vals_c2 = (double *)malloc(nc * sizeof(double)); - double *errors_c2 = (double *)malloc(nc * sizeof(double)); - double *vals_c3 = (double *)malloc(nc * sizeof(double)); - double *errors_c3 = (double *)malloc(nc * sizeof(double)); - double *vals_c4 = (double *)malloc(nc * sizeof(double)); - double *errors_c4 = (double *)malloc(nc * sizeof(double)); - - char *out_filename = (char *)malloc(100 * sizeof(char)); - uint_t out_filename_len = 0; - - for (uint_t i = 0; i < nc; i++) { - char *fn = argv[1 + i]; - char *buff = (char *)malloc(20 * sizeof(char)); - uint_t pos = 0; - uint_t c = 0; - while (c < 5) { - if (fn[pos] == '_') { - c++; - } - if (i == 0) { - out_filename[pos] = fn[pos]; - out_filename_len++; - } - pos++; - } - c = 0; - while (fn[pos] != '_') { - buff[c] = fn[pos]; - pos++; - c++; - } - buff[c] = '\0'; - Ls_c[i] = atoi(buff); - c = 0; - pos++; - while (fn[pos] != '_') { - buff[c] = fn[pos]; - pos++; - c++; - } - buff[c] = '\0'; - betas_c[i] = atof(buff); - free(buff); - - uint_t dist_len = 4 * pow(Ls_c[i], 2); - uint32_t *dist = malloc(dist_len * sizeof(uint32_t)); - FILE *dist_file = fopen(fn, "rb"); - fread(dist, sizeof(uint32_t), dist_len, dist_file); - fclose(dist_file); - { - double mom1[2]; - mom(dist_len, 1, dist, mom1); - vals_c1[i] = mom1[0]; - errors_c1[i] = mom1[1]; - - double mom2[2]; - mom(dist_len, 2, dist, mom2); - vals_c2[i] = mom2[0]; - errors_c2[i] = mom2[1]; - - double mom3[2]; - mom(dist_len, 3, dist, mom3); - vals_c3[i] = mom3[0]; - errors_c3[i] = mom3[1]; - - double mom4[2]; - mom(dist_len, 4, dist, mom4); - vals_c4[i] = mom4[0]; - errors_c4[i] = mom4[1]; - } - free(dist); - } - - out_filename[out_filename_len - 1] = '.'; - out_filename[out_filename_len] = 't'; - out_filename[out_filename_len + 1] = 'x'; - out_filename[out_filename_len + 2] = 't'; - out_filename[out_filename_len + 3] = '\0'; - - FILE *out_file = fopen(out_filename, "w"); - - for (uint_t i = 0; i < nc; i++) { - fprintf(out_file, "%u %g %g %g %g %g %g %g %g %g\n", Ls_c[i], betas_c[i], - vals_c1[i], errors_c1[i], vals_c2[i], errors_c2[i], vals_c3[i], - errors_c3[i], vals_c4[i], errors_c4[i]); - } - - fclose(out_file); - - free(Ls_c); - free(betas_c); - free(vals_c1); - free(errors_c1); - free(vals_c2); - free(errors_c2); - free(vals_c3); - free(errors_c3); - - return 0; -} diff --git a/src/big_anal_process.c b/src/big_anal_process.c deleted file mode 100644 index 8c1d8ba..0000000 --- a/src/big_anal_process.c +++ /dev/null @@ -1,158 +0,0 @@ - -#include "fracture.h" -#include <gsl/gsl_blas.h> -#include <gsl/gsl_matrix.h> -#include <gsl/gsl_sf_erf.h> -#include <gsl/gsl_sf_laguerre.h> -#include <gsl/gsl_vector.h> -#include <sys/stat.h> - -void get_mean(uint_t len, double *data, double result[2]) { - double total = 0; - - for (uint_t i = 0; i < len; i++) { - total += data[i]; - } - - double mean = total / len; - double total_sq_diff = 0; - - for (uint_t i = 0; i < len; i++) { - total_sq_diff += pow(mean - data[i], 2); - } - - double mean_err = sqrt(total_sq_diff) / len; - - result[0] = mean; - result[1] = mean_err; -} - -void get_mom(uint_t len, uint_t n, double *data, double mean[2], - double result[2]) { - double total_n_diff = 0; - double total_n2_diff = 0; - - for (uint_t i = 0; i < len; i++) { - total_n_diff += pow(fabs(mean[0] - data[i]), n); - total_n2_diff += pow(fabs(mean[0] - data[i]), 2 * n); - } - - double mom = total_n_diff / len; - double mom_err = sqrt(total_n2_diff) / len; - - result[0] = mom; - result[1] = mom_err; -} - -int main(int argc, char *argv[]) { - uint_t nc = argc - 1; - uint_t *Ls_c = (uint_t *)malloc(nc * sizeof(uint_t)); - double *betas_c = (double *)malloc(nc * sizeof(double)); - double *vals_c1 = (double *)malloc(nc * sizeof(double)); - double *errors_c1 = (double *)malloc(nc * sizeof(double)); - double *vals_c2 = (double *)malloc(nc * sizeof(double)); - double *errors_c2 = (double *)malloc(nc * sizeof(double)); - double *vals_c3 = (double *)malloc(nc * sizeof(double)); - double *errors_c3 = (double *)malloc(nc * sizeof(double)); - double *vals_c4 = (double *)malloc(nc * sizeof(double)); - double *errors_c4 = (double *)malloc(nc * sizeof(double)); - - char *out_filename = (char *)malloc(100 * sizeof(char)); - uint_t out_filename_len = 0; - - for (uint_t i = 0; i < nc; i++) { - char *fn = argv[1 + i]; - char *buff = (char *)malloc(20 * sizeof(char)); - uint_t pos = 0; - uint_t c = 0; - while (c < 5) { - if (fn[pos] == '_') { - c++; - } - if (i == 0) { - out_filename[pos] = fn[pos]; - out_filename_len++; - } - pos++; - } - c = 0; - while (fn[pos] != '_') { - buff[c] = fn[pos]; - pos++; - c++; - } - buff[c] = '\0'; - Ls_c[i] = atoi(buff); - c = 0; - pos++; - while (fn[pos] != '_') { - buff[c] = fn[pos]; - pos++; - c++; - } - buff[c] = '\0'; - betas_c[i] = atof(buff); - free(buff); - - struct stat info; - stat(fn, &info); - uint_t num_bytes = info.st_size; - uint_t dist_len = (num_bytes * sizeof(char)) / sizeof(double); - - double *dist = malloc(dist_len * sizeof(double)); - FILE *dist_file = fopen(fn, "rb"); - fread(dist, sizeof(double), dist_len, dist_file); - fclose(dist_file); - { - double mom1[2]; - get_mean(dist_len, dist, mom1); - vals_c1[i] = mom1[0]; - errors_c1[i] = mom1[1]; - - double mom2[2]; - get_mom(dist_len, 2, dist, mom1, mom2); - vals_c2[i] = mom2[0]; - errors_c2[i] = mom2[1]; - - double mom3[2]; - get_mom(dist_len, 3, dist, mom1, mom3); - vals_c3[i] = mom3[0]; - errors_c3[i] = mom3[1]; - - double mom4[2]; - get_mom(dist_len, 4, dist, mom1, mom4); - vals_c4[i] = mom4[0]; - errors_c4[i] = mom4[1]; - } - free(dist); - } - - out_filename[out_filename_len - 1] = '.'; - out_filename[out_filename_len] = 't'; - out_filename[out_filename_len + 1] = 'x'; - out_filename[out_filename_len + 2] = 't'; - out_filename[out_filename_len + 3] = '\0'; - - FILE *out_file = fopen(out_filename, "w"); - - for (uint_t i = 0; i < nc; i++) { - fprintf(out_file, "%u %g %g %g %g %g %g %g %g %g\n", Ls_c[i], betas_c[i], - vals_c1[i], errors_c1[i], vals_c2[i], errors_c2[i], vals_c3[i], - errors_c3[i], vals_c4[i], errors_c4[i]); - } - - fclose(out_file); - - free(Ls_c); - free(betas_c); - free(vals_c1); - free(errors_c1); - free(vals_c2); - free(errors_c2); - free(vals_c3); - free(errors_c3); - free(vals_c4); - free(errors_c4); - - return 0; -} diff --git a/src/cen_anal_process.c b/src/cen_anal_process.c deleted file mode 100644 index ee2b029..0000000 --- a/src/cen_anal_process.c +++ /dev/null @@ -1,157 +0,0 @@ - -#include "fracture.h" -#include <gsl/gsl_blas.h> -#include <gsl/gsl_matrix.h> -#include <gsl/gsl_sf_erf.h> -#include <gsl/gsl_sf_laguerre.h> -#include <gsl/gsl_vector.h> - -void get_mean(uint_t len, uint32_t *data, double result[2]) { - uint_t total = 0; - double mean = 0; - double mean_err = 0; - - for (uint_t i = 0; i < len; i++) { - uint32_t datai = data[i]; - - total += datai; - mean += i * datai; - mean_err += gsl_pow_2(i) * datai; - } - - double meanf = mean / total; - double meanf_err = meanf * sqrt(mean_err / gsl_pow_2(mean) + 1 / total); - - result[0] = meanf; - result[1] = meanf_err; -} - -void get_mom(uint_t len, uint_t n, uint32_t *data, double mean[2], - double result[2]) { - uint_t total = 0; - double mom = 0; - double mom_err = 0; - - for (uint_t i = 0; i < len; i++) { - uint32_t datai = data[i]; - double in = pow(fabs(((double)i) - mean[0]), n); - - total += datai; - mom += in * datai; - mom_err += gsl_pow_2(in) * - datai; // + gsl_pow_2(n * mean[1] / (((double)i) - mean[0]))); - } - - double momf = mom / total; - double momf_err = momf * sqrt(mom_err / gsl_pow_2(mom) + 1 / total); - - result[0] = momf; - result[1] = momf_err; -} - -int main(int argc, char *argv[]) { - uint_t nc = argc - 1; - uint_t *Ls_c = (uint_t *)malloc(nc * sizeof(uint_t)); - double *betas_c = (double *)malloc(nc * sizeof(double)); - double *vals_c1 = (double *)malloc(nc * sizeof(double)); - double *errors_c1 = (double *)malloc(nc * sizeof(double)); - double *vals_c2 = (double *)malloc(nc * sizeof(double)); - double *errors_c2 = (double *)malloc(nc * sizeof(double)); - double *vals_c3 = (double *)malloc(nc * sizeof(double)); - double *errors_c3 = (double *)malloc(nc * sizeof(double)); - double *vals_c4 = (double *)malloc(nc * sizeof(double)); - double *errors_c4 = (double *)malloc(nc * sizeof(double)); - - char *out_filename = (char *)malloc(100 * sizeof(char)); - uint_t out_filename_len = 0; - - for (uint_t i = 0; i < nc; i++) { - char *fn = argv[1 + i]; - char *buff = (char *)malloc(20 * sizeof(char)); - uint_t pos = 0; - uint_t c = 0; - while (c < 5) { - if (fn[pos] == '_') { - c++; - } - if (i == 0) { - out_filename[pos] = fn[pos]; - out_filename_len++; - } - pos++; - } - c = 0; - while (fn[pos] != '_') { - buff[c] = fn[pos]; - pos++; - c++; - } - buff[c] = '\0'; - Ls_c[i] = atoi(buff); - c = 0; - pos++; - while (fn[pos] != '_') { - buff[c] = fn[pos]; - pos++; - c++; - } - buff[c] = '\0'; - betas_c[i] = atof(buff); - free(buff); - - uint_t dist_len = 4 * pow(Ls_c[i], 2); - uint32_t *dist = malloc(dist_len * sizeof(uint32_t)); - FILE *dist_file = fopen(fn, "rb"); - fread(dist, sizeof(uint32_t), dist_len, dist_file); - fclose(dist_file); - { - double mom1[2]; - get_mean(dist_len, dist, mom1); - vals_c1[i] = mom1[0]; - errors_c1[i] = mom1[1]; - - double mom2[2]; - get_mom(dist_len, 2, dist, mom1, mom2); - vals_c2[i] = mom2[0]; - errors_c2[i] = mom2[1]; - - double mom3[2]; - get_mom(dist_len, 3, dist, mom1, mom3); - vals_c3[i] = mom3[0]; - errors_c3[i] = mom3[1]; - - double mom4[2]; - get_mom(dist_len, 4, dist, mom1, mom4); - vals_c4[i] = mom4[0]; - errors_c4[i] = mom4[1]; - } - free(dist); - } - - out_filename[out_filename_len - 1] = '.'; - out_filename[out_filename_len] = 't'; - out_filename[out_filename_len + 1] = 'x'; - out_filename[out_filename_len + 2] = 't'; - out_filename[out_filename_len + 3] = '\0'; - - FILE *out_file = fopen(out_filename, "w"); - - for (uint_t i = 0; i < nc; i++) { - fprintf(out_file, "%u %g %g %g %g %g %g %g %g %g\n", Ls_c[i], betas_c[i], - vals_c1[i], errors_c1[i], vals_c2[i], errors_c2[i], vals_c3[i], - errors_c3[i], vals_c4[i], errors_c4[i]); - } - - fclose(out_file); - - free(Ls_c); - free(betas_c); - free(vals_c1); - free(errors_c1); - free(vals_c2); - free(errors_c2); - free(vals_c3); - free(errors_c3); - - return 0; -} diff --git a/src/corr_test.c b/src/corr_test.c deleted file mode 100644 index cb00212..0000000 --- a/src/corr_test.c +++ /dev/null @@ -1,27 +0,0 @@ - -#include "fracture.h" - -int main() { - cholmod_common c; - CHOL_F(start)(&c); - - unsigned int width = 64; - - graph_t *network = graph_create(VORONOI_LATTICE, TORUS_BOUND, 128, false); - net_t *instance = net_create(network, 1e-14, 3, 0, true, &c); - net_fracture(instance, &c, 1e-10); - - double *corr = get_corr(instance, NULL, &c); - - for (int i = 0; i < 2 * width; i++) { - printf("%g ", corr[i]); - } - printf("\n"); - - net_free(instance, &c); - graph_free(network); - - CHOL_F(finish)(&c); - - return 0; -} diff --git a/src/fracture.c b/src/fracture.c deleted file mode 100644 index 5b44238..0000000 --- a/src/fracture.c +++ /dev/null @@ -1,1068 +0,0 @@ - -#include <fftw3.h> - -#include "fracture.h" - -int main(int argc, char *argv[]) { - int opt; - - fftw_set_timelimit(1); - - // defining variables to be (potentially) set by command line flags - uint8_t filename_len; - uint32_t N; - uint_t L; - double beta, inf, cutoff, crack_len; - bool save_data, save_cluster_dist, use_voltage_boundaries, use_dual, - save_network, save_crit_stress, save_energy, save_conductivity, - save_damage, save_threshold, save_current_load, save_correlations; - bound_t boundary; - lattice_t lattice; - - // assume filenames will be less than 100 characters - - filename_len = 100; - - // set default values - - N = 100; - L = 16; - crack_len = 0.; - beta = .3; - inf = 1e10; - cutoff = 1e-9; - boundary = FREE_BOUND; - lattice = VORONOI_LATTICE; - save_data = false; - save_cluster_dist = false; - use_voltage_boundaries = false; - use_dual = false; - save_network = false; - save_crit_stress = false; - save_damage = false; - save_conductivity = false; - save_energy = false; - save_threshold = false; - save_current_load = false; - save_correlations = false; - - uint8_t bound_i; - char boundc2 = 'f'; - uint8_t lattice_i; - char lattice_c = 'v'; - char dual_c = 'o'; - - // get commandline options - - while ((opt = getopt(argc, argv, "n:L:b:B:q:dVcoNsCrDl:TEz")) != -1) { - switch (opt) { - case 'n': - N = atoi(optarg); - break; - case 'L': - L = atoi(optarg); - break; - case 'b': - beta = atof(optarg); - break; - case 'l': - crack_len = atof(optarg); - break; - case 'B': - bound_i = atoi(optarg); - switch (bound_i) { - case 0: - boundary = FREE_BOUND; - boundc2 = 'f'; - break; - case 1: - boundary = CYLINDER_BOUND; - boundc2 = 'c'; - break; - case 2: - boundary = TORUS_BOUND; - use_voltage_boundaries = true; - boundc2 = 't'; - break; - case 3: - boundary = EMBEDDED_BOUND; - boundc2 = 'e'; - use_dual = true; - use_voltage_boundaries = true; - break; - default: - printf("boundary specifier must be 0 (FREE_BOUND), 1 (CYLINDER_BOUND), " - "or 2 (TORUS_BOUND).\n"); - exit(EXIT_FAILURE); - } - break; - case 'q': - lattice_i = atoi(optarg); - switch (lattice_i) { - case 0: - lattice = VORONOI_LATTICE; - lattice_c = 'v'; - break; - case 1: - lattice = DIAGONAL_LATTICE; - lattice_c = 'd'; - break; - case 2: - lattice = VORONOI_HYPERUNIFORM_LATTICE; - lattice_c = 'h'; - break; - case 3: - lattice = TRIANGULAR_LATTICE; - lattice_c = 't'; - break; - case 4: - lattice = SQUARE_LATTICE; - lattice_c = 's'; - break; - default: - printf("lattice specifier must be 0 (VORONOI_LATTICE), 1 " - "(DIAGONAL_LATTICE), or 2 (VORONOI_HYPERUNIFORM_LATTICE).\n"); - exit(EXIT_FAILURE); - } - break; - case 'd': - save_damage = true; - break; - case 'V': - use_voltage_boundaries = true; - break; - case 'D': - use_dual = true; - dual_c = 'd'; - break; - case 'c': - save_cluster_dist = true; - break; - case 'o': - save_data = true; - break; - case 'N': - save_network = true; - break; - case 's': - save_crit_stress = true; - break; - case 'r': - save_conductivity = true; - break; - case 'E': - save_energy = true; - break; - case 'T': - save_threshold = true; - break; - case 'C': - save_current_load = true; - break; - case 'z': - save_correlations = true; - break; - default: /* '?' */ - exit(EXIT_FAILURE); - } - } - - char boundc; - if (use_voltage_boundaries) { - boundc = 'v'; - } else { - boundc = 'c'; - } - - double *dd_correlations; // damage-damage - double *dc_correlations; // damage-crack - double *db_correlations; // damage-backbone - double *ds_correlations; // damage-stress - double *dA_correlations; // damage-avalanche - double *cc_correlations; // crack-crack - double *cb_correlations; // crack-backbone - double *cs_correlations; // crack-stress - double *cA_correlations; // crack-avalanche - double *bb_correlations; // backbone-backbone - double *bs_correlations; // backbone-stress - double *bA_correlations; // backbone-avalanche - double *ss_correlations; // stress-stress - double *AA_correlations; // avalanche-avalanche - double *DD_correlations; // after-crack damage-damage - double *DC_correlations; // after-crack damage-crack - double *DB_correlations; // after-crack damage-backbone - double *CC_correlations; // after-crack crack-crack - double *CB_correlations; // after-crack crack-backbone - double *BB_correlations; // after-crack backbone-backbone - double *fftw_forward_in; - fftw_complex *fftw_forward_out; - fftw_complex *fftw_reverse_in; - double *fftw_reverse_out; - fftw_plan forward_plan; - fftw_plan reverse_plan; - uint64_t N_averaged = 0; - double mean_D = 0; - double mean_A = 0; - double mean_B = 0; - double mean_C = 0; - double mean_d = 0; - double mean_c = 0; - double mean_b = 0; - double mean_s = 0; - char *correlations_filename; - if (save_correlations) { - assert(lattice == DIAGONAL_LATTICE); - dd_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - dc_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - db_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - ds_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - dA_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - cc_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - cb_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - cs_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - cA_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - bb_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - bs_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - bA_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - ss_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - AA_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - DD_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - DC_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - DB_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - CC_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - CB_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - BB_correlations = (double *)calloc(pow(L, 2), sizeof(double)); - fftw_forward_in = (double *)fftw_malloc(pow(L, 2) * sizeof(double)); - fftw_forward_out = (fftw_complex *)fftw_malloc(pow(L, 2) * sizeof(fftw_complex)); - fftw_reverse_in = (fftw_complex *)fftw_malloc(pow(L, 2) * sizeof(fftw_complex)); - fftw_reverse_out = (double *)fftw_malloc(pow(L, 2) * sizeof(double)); - forward_plan = fftw_plan_dft_r2c_2d(L, L, fftw_forward_in, fftw_forward_out, 0); - reverse_plan = fftw_plan_dft_c2r_2d(L, L, fftw_reverse_in, fftw_reverse_out, 0); - - correlations_filename = (char *)malloc(filename_len * sizeof(char)); - snprintf(correlations_filename, filename_len, "corr_%c_%c_%c_%c_%d_%g_%g.dat", - lattice_c, dual_c, boundc, boundc2, L, beta, crack_len); - - FILE *correlations_out = fopen(correlations_filename, "rb"); - - if (correlations_out != NULL) { - fread(&N_averaged, sizeof(uint64_t), 1, correlations_out); - fread(&mean_d, sizeof(double), 1, correlations_out); - fread(&mean_c, sizeof(double), 1, correlations_out); - fread(&mean_b, sizeof(double), 1, correlations_out); - fread(&mean_s, sizeof(double), 1, correlations_out); - fread(&mean_A, sizeof(double), 1, correlations_out); - fread(&mean_D, sizeof(double), 1, correlations_out); - fread(&mean_C, sizeof(double), 1, correlations_out); - fread(&mean_B, sizeof(double), 1, correlations_out); - fread(dd_correlations, sizeof(double), pow(L, 2), correlations_out); - fread(dc_correlations, sizeof(double), pow(L, 2), correlations_out); - fread(db_correlations, sizeof(double), pow(L, 2), correlations_out); - fread(ds_correlations, sizeof(double), pow(L, 2), correlations_out); - fread(dA_correlations, sizeof(double), pow(L, 2), correlations_out); - fread(cc_correlations, sizeof(double), pow(L, 2), correlations_out); - fread(cb_correlations, sizeof(double), pow(L, 2), correlations_out); - fread(cs_correlations, sizeof(double), pow(L, 2), correlations_out); - fread(cA_correlations, sizeof(double), pow(L, 2), correlations_out); - fread(bb_correlations, sizeof(double), pow(L, 2), correlations_out); - fread(bs_correlations, sizeof(double), pow(L, 2), correlations_out); - fread(bA_correlations, sizeof(double), pow(L, 2), correlations_out); - fread(ss_correlations, sizeof(double), pow(L, 2), correlations_out); - fread(AA_correlations, sizeof(double), pow(L, 2), correlations_out); - fread(DD_correlations, sizeof(double), pow(L, 2), correlations_out); - fread(DC_correlations, sizeof(double), pow(L, 2), correlations_out); - fread(DB_correlations, sizeof(double), pow(L, 2), correlations_out); - fread(CC_correlations, sizeof(double), pow(L, 2), correlations_out); - fread(CB_correlations, sizeof(double), pow(L, 2), correlations_out); - fread(BB_correlations, sizeof(double), pow(L, 2), correlations_out); - fclose(correlations_out); - } - } - - FILE *data_out; - - if (save_data) { - char *data_filename = (char *)malloc(filename_len * sizeof(char)); - snprintf(data_filename, filename_len, "data_%c_%c_%c_%c_%u_%g_%g.txt", - lattice_c, dual_c, boundc, boundc2, L, beta, crack_len); - data_out = fopen(data_filename, "a"); - free(data_filename); - } - - uint_t max_verts, max_edges; - - // these are very liberal estimates - max_verts = 4 * pow(L, 2); - max_edges = 4 * pow(L, 2); - - if (max_verts > CINT_MAX) { - exit(EXIT_FAILURE); - } - - // define arrays for saving cluster and avalanche distributions - uint32_t *cluster_size_dist; - uint32_t *avalanche_size_dist; - char *c_filename; - char *a_filename; - if (save_cluster_dist) { - cluster_size_dist = (uint32_t *)calloc(max_verts, sizeof(uint32_t)); - avalanche_size_dist = (uint32_t *)calloc(max_edges, sizeof(uint32_t)); - - c_filename = (char *)malloc(filename_len * sizeof(char)); - a_filename = (char *)malloc(filename_len * sizeof(char)); - snprintf(c_filename, filename_len, "cstr_%c_%c_%c_%c_%d_%g_%g.dat", - lattice_c, dual_c, boundc, boundc2, L, beta, crack_len); - snprintf(a_filename, filename_len, "avln_%c_%c_%c_%c_%d_%g_%g.dat", - lattice_c, dual_c, boundc, boundc2, L, beta, crack_len); - - FILE *cluster_out = fopen(c_filename, "rb"); - FILE *avalanche_out = fopen(a_filename, "rb"); - - if (cluster_out != NULL) { - fread(cluster_size_dist, sizeof(uint32_t), max_verts, cluster_out); - fclose(cluster_out); - } - if (avalanche_out != NULL) { - fread(avalanche_size_dist, sizeof(uint32_t), max_edges, avalanche_out); - fclose(avalanche_out); - } - } - - long double *crit_stress; - if (save_crit_stress) { - crit_stress = (long double *)malloc(N * sizeof(long double)); - } - - double *conductivity; - if (save_conductivity) { - conductivity = (double *)malloc(N * sizeof(double)); - } - - // define arrays for saving damage distributions - uint32_t *damage; - char *d_filename; - if (save_damage) { - damage = (uint32_t *)calloc(max_edges, sizeof(uint32_t)); - - d_filename = (char *)malloc(filename_len * sizeof(char)); - snprintf(d_filename, filename_len, "damg_%c_%c_%c_%c_%d_%g_%g.dat", - lattice_c, dual_c, boundc, boundc2, L, beta, crack_len); - - FILE *damage_out = fopen(d_filename, "rb"); - - if (damage_out != NULL) { - fread(damage, sizeof(uint32_t), max_edges, damage_out); - fclose(damage_out); - } - } - - long double *energy; - if (save_energy) { - energy = (long double *)malloc(N * sizeof(long double)); - } - - long double *thresholds; - if (save_threshold) { - thresholds = (long double *)malloc(N * sizeof(long double)); - } - - long double *loads; - if (save_current_load) { - loads = (long double *)malloc(N * sizeof(long double)); - } - - // start cholmod - cholmod_common c; - CHOL_F(start)(&c); - - /* if we use voltage boundary conditions, the laplacian matrix is positive - * definite and we can use a supernodal LL decomposition. otherwise we need - * to use the simplicial LDL decomposition - */ - if (use_voltage_boundaries) { - //(&c)->supernodal = CHOLMOD_SUPERNODAL; - (&c)->supernodal = CHOLMOD_SIMPLICIAL; - } else { - (&c)->supernodal = CHOLMOD_SIMPLICIAL; - } - - printf("\n"); - for (uint32_t i = 0; i < N; i++) { - printf("\033[F\033[JFRACTURE: %0*d / %d\n", (uint8_t)log10(N) + 1, i + 1, - N); - - graph_t *g = graph_create(lattice, boundary, L, use_dual); - net_t *net = - net_create(g, inf, beta, crack_len, use_voltage_boundaries, &c); - net_t *tmp_net = net_copy(net, &c); - data_t *data = net_fracture(tmp_net, &c, cutoff); - - uint_t max_pos = 0; - long double max_val = 0; - - double cond0; - { - double *tmp_voltages = net_voltages(net, &c); - cond0 = net_conductivity(net, tmp_voltages); - free(tmp_voltages); - } - - for (uint_t j = 0; j < data->num_broken; j++) { - long double val = data->extern_field[j]; - - if (val > max_val) { - max_pos = j; - max_val = val; - } - } - - uint_t av_size = 0; - long double cur_val = 0; - - for (uint_t j = 0; j < max_pos; j++) { - uint_t next_broken = data->break_list[j]; - - break_edge(net, next_broken, &c); - - long double val = data->extern_field[j]; - if (save_cluster_dist) { - if (val < cur_val) { - av_size++; - } - - if (val > cur_val) { - avalanche_size_dist[av_size]++; - av_size = 0; - cur_val = val; - } - } - } - - if (save_correlations) { - uint32_t damage1 = 0; - for (uint32_t j = 0; j < g->ne; j++) { - if (tmp_net->fuses[j]) { - fftw_forward_in[j] = 1.0; - damage1++; - } else { - fftw_forward_in[j] = 0.0; - } - } - - fftw_execute(forward_plan); - fftw_complex *D_transform = (fftw_complex *)malloc(pow(L, 2) * sizeof(fftw_complex)); - memcpy(D_transform, fftw_forward_out, g->ne * sizeof(fftw_complex)); - - dll_t *cycle = find_cycles(g, tmp_net->fuses); - components_t *comp = get_clusters(tmp_net); - - uint32_t in_crack1 = 0; - for (uint32_t j = 0; j < g->ne; j++) { - if (tmp_net->fuses[j] && comp->labels[g->dev[2 * cycle->e]] == comp->labels[g->dev[2 * j]]) { - fftw_forward_in[j] = 1.0; - in_crack1++; - } else { - fftw_forward_in[j] = 0.0; - } - } - - fftw_execute(forward_plan); - fftw_complex *C_transform = (fftw_complex *)malloc(pow(L, 2) * sizeof(fftw_complex)); - memcpy(C_transform, fftw_forward_out, g->ne * sizeof(fftw_complex)); - - for (uint32_t j = 0; j < g->ne; j++) { - fftw_forward_in[j] = 0.0; - } - - uint32_t in_backbone1 = 0; - dll_t *tmp_cycle = cycle; - while (tmp_cycle != NULL) { - fftw_forward_in[tmp_cycle->e] = 1.0; - in_backbone1++; - tmp_cycle = tmp_cycle->right; - } - - fftw_execute(forward_plan); - fftw_complex *B_transform = (fftw_complex *)malloc(g->ne * sizeof(fftw_complex)); - memcpy(B_transform, fftw_forward_out, g->ne * sizeof(fftw_complex)); - - uint32_t in_avalanche = 0; - for (uint32_t j = 0; j < g->ne; j++) { - if (tmp_net->fuses[j] && !(net->fuses[j])) { - fftw_forward_in[j] = 1.0; - } else { - fftw_forward_in[j] = 0.0; - } - } - - fftw_execute(forward_plan); - fftw_complex *A_transform = (fftw_complex *)malloc(g->ne * sizeof(fftw_complex)); - memcpy(A_transform, fftw_forward_out, g->ne * sizeof(fftw_complex)); - - uint32_t damage2 = 0; - for (uint32_t j = 0; j < g->ne; j++) { - if (net->fuses[j]) { - fftw_forward_in[j] = 1.0; - damage2++; - } else { - fftw_forward_in[j] = 0.0; - } - } - - fftw_execute(forward_plan); - fftw_complex *d_transform = (fftw_complex *)malloc(pow(L, 2) * sizeof(fftw_complex)); - memcpy(d_transform, fftw_forward_out, g->ne * sizeof(fftw_complex)); - - uint32_t in_crack2 = 0; - for (uint32_t j = 0; j < g->ne; j++) { - if (net->fuses[j] && comp->labels[2 * g->dev[cycle->e]] == comp->labels[g->dev[2 * j]]) { - fftw_forward_in[j] = 1.0; - in_crack2++; - } else { - fftw_forward_in[j] = 0.0; - } - } - - fftw_execute(forward_plan); - fftw_complex *c_transform = (fftw_complex *)malloc(g->ne * sizeof(fftw_complex)); - memcpy(c_transform, fftw_forward_out, g->ne * sizeof(fftw_complex)); - - for (uint32_t j = 0; j < g->ne; j++) { - fftw_forward_in[j] = 0.0; - } - - uint32_t in_backbone2 = 0; - tmp_cycle = cycle; - while (tmp_cycle != NULL) { - if (net->fuses[tmp_cycle->e]) { - fftw_forward_in[tmp_cycle->e] = 1.0; - in_backbone2++; - } - tmp_cycle = tmp_cycle->right; - } - - fftw_execute(forward_plan); - fftw_complex *b_transform = (fftw_complex *)malloc(g->ne * sizeof(fftw_complex)); - memcpy(b_transform, fftw_forward_out, g->ne * sizeof(fftw_complex)); - - graph_components_free(comp); - while (cycle != NULL) { - dll_t *old = cycle; - cycle = cycle->right; - free(old); - } - - double *tmp_voltage = net_voltages(net, &c); - double *tmp_current = net_currents(net, tmp_voltage, &c); - free(tmp_voltage); - - double t_stress = 0; - for (uint32_t j = 0; j < g->ne; j++) { - fftw_forward_in[j] = tmp_current[j]; - t_stress += tmp_current[j]; - } - - free(tmp_current); - - fftw_execute(forward_plan); - fftw_complex *s_transform = (fftw_complex *)malloc(pow(L, 2) * sizeof(fftw_complex)); - memcpy(s_transform, fftw_forward_out, g->ne * sizeof(fftw_complex)); - - mean_D = (double)damage1 / (1.0 + N_averaged) + (double)N_averaged * mean_D / (N_averaged + 1.0); - mean_C = (double)in_crack1 / (1.0 + N_averaged) + (double)N_averaged * mean_C / (N_averaged + 1.0); - mean_B = (double)in_backbone1 / (1.0 + N_averaged) + (double)N_averaged * mean_B / (N_averaged + 1.0); - mean_A = (double)in_avalanche / (1.0 + N_averaged) + (double)N_averaged * mean_A / (N_averaged + 1.0); - mean_d = (double)damage2 / (1.0 + N_averaged) + (double)N_averaged * mean_d / (N_averaged + 1.0); - mean_c = (double)in_crack2 / (1.0 + N_averaged) + (double)N_averaged * mean_c / (N_averaged + 1.0); - mean_b = (double)in_backbone2 / (1.0 + N_averaged) + (double)N_averaged * mean_b / (N_averaged + 1.0); - mean_s = (double)t_stress / (1.0 + N_averaged) + (double)N_averaged * mean_s / (N_averaged + 1.0); - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = D_transform[j][0] * D_transform[j][0] + D_transform[j][1] * D_transform[j][1]; - fftw_reverse_in[j][1] = 0.0; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - DD_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * DD_correlations[j] / (N_averaged + 1.0); - } - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = D_transform[j][0] * C_transform[j][0] + D_transform[j][1] * C_transform[j][1]; - fftw_reverse_in[j][1] = D_transform[j][0] * C_transform[j][1] - D_transform[j][1] * C_transform[j][0]; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - DC_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * DC_correlations[j] / (N_averaged + 1.0); - } - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = D_transform[j][0] * B_transform[j][0] + D_transform[j][1] * B_transform[j][1]; - fftw_reverse_in[j][1] = D_transform[j][0] * B_transform[j][1] - D_transform[j][1] * B_transform[j][0]; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - DB_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * DB_correlations[j] / (N_averaged + 1.0); - } - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = C_transform[j][0] * C_transform[j][0] + C_transform[j][1] * C_transform[j][1]; - fftw_reverse_in[j][1] = 0.0; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - CC_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * CC_correlations[j] / (N_averaged + 1.0); - } - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = C_transform[j][0] * B_transform[j][0] + C_transform[j][1] * B_transform[j][1]; - fftw_reverse_in[j][1] = C_transform[j][0] * B_transform[j][1] - C_transform[j][1] * B_transform[j][0]; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - CB_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * CB_correlations[j] / (N_averaged + 1.0); - } - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = B_transform[j][0] * B_transform[j][0] + B_transform[j][1] * B_transform[j][1]; - fftw_reverse_in[j][1] = 0.0; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - BB_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * BB_correlations[j] / (N_averaged + 1.0); - } - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = A_transform[j][0] * A_transform[j][0] + A_transform[j][1] * A_transform[j][1]; - fftw_reverse_in[j][1] = 0.0; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - AA_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * AA_correlations[j] / (N_averaged + 1.0); - } - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = d_transform[j][0] * d_transform[j][0] + d_transform[j][1] * d_transform[j][1]; - fftw_reverse_in[j][1] = 0.0; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - dd_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * dd_correlations[j] / (N_averaged + 1.0); - } - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = d_transform[j][0] * c_transform[j][0] + d_transform[j][1] * c_transform[j][1]; - fftw_reverse_in[j][1] = d_transform[j][0] * c_transform[j][1] - d_transform[j][1] * c_transform[j][0]; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - dc_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * dc_correlations[j] / (N_averaged + 1.0); - } - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = d_transform[j][0] * s_transform[j][0] + d_transform[j][1] * s_transform[j][1]; - fftw_reverse_in[j][1] = d_transform[j][0] * s_transform[j][1] - d_transform[j][1] * s_transform[j][0]; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - ds_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * ds_correlations[j] / (N_averaged + 1.0); - } - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = d_transform[j][0] * b_transform[j][0] + d_transform[j][1] * b_transform[j][1]; - fftw_reverse_in[j][1] = d_transform[j][0] * b_transform[j][1] - d_transform[j][1] * b_transform[j][0]; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - db_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * db_correlations[j] / (N_averaged + 1.0); - } - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = d_transform[j][0] * A_transform[j][0] + d_transform[j][1] * A_transform[j][1]; - fftw_reverse_in[j][1] = d_transform[j][0] * A_transform[j][1] - d_transform[j][1] * A_transform[j][0]; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - dA_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * dA_correlations[j] / (N_averaged + 1.0); - } - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = c_transform[j][0] * c_transform[j][0] + c_transform[j][1] * c_transform[j][1]; - fftw_reverse_in[j][1] = 0.0; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - cc_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * cc_correlations[j] / (N_averaged + 1.0); - } - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = c_transform[j][0] * s_transform[j][0] + c_transform[j][1] * s_transform[j][1]; - fftw_reverse_in[j][1] = c_transform[j][0] * s_transform[j][1] - c_transform[j][1] * s_transform[j][0]; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - cs_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * cs_correlations[j] / (N_averaged + 1.0); - } - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = c_transform[j][0] * b_transform[j][0] + c_transform[j][1] * b_transform[j][1]; - fftw_reverse_in[j][1] = c_transform[j][0] * b_transform[j][1] - c_transform[j][1] * b_transform[j][0]; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - cb_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * cb_correlations[j] / (N_averaged + 1.0); - } - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = c_transform[j][0] * A_transform[j][0] + c_transform[j][1] * A_transform[j][1]; - fftw_reverse_in[j][1] = c_transform[j][0] * A_transform[j][1] - c_transform[j][1] * A_transform[j][0]; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - cA_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * cA_correlations[j] / (N_averaged + 1.0); - } - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = b_transform[j][0] * b_transform[j][0] + b_transform[j][1] * b_transform[j][1]; - fftw_reverse_in[j][1] = 0.0; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - bb_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * bb_correlations[j] / (N_averaged + 1.0); - } - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = b_transform[j][0] * s_transform[j][0] + s_transform[j][1] * s_transform[j][1]; - fftw_reverse_in[j][1] = b_transform[j][0] * s_transform[j][1] - s_transform[j][1] * s_transform[j][0]; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - bs_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * bs_correlations[j] / (N_averaged + 1.0); - } - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = A_transform[j][0] * b_transform[j][0] + A_transform[j][1] * b_transform[j][1]; - fftw_reverse_in[j][1] = A_transform[j][0] * b_transform[j][1] - A_transform[j][1] * b_transform[j][0]; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - bA_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * bA_correlations[j] / (N_averaged + 1.0); - } - - - for (uint32_t j = 0; j < L * (L / 2 + 1); j++) { - fftw_reverse_in[j][0] = s_transform[j][0] * s_transform[j][0] + s_transform[j][1] * s_transform[j][1]; - fftw_reverse_in[j][1] = 0.0; - } - - fftw_execute(reverse_plan); - - for (uint32_t j = 0; j < g->ne; j++) { - ss_correlations[j] = fftw_reverse_out[j] / (1.0 + N_averaged) + (double)N_averaged * ss_correlations[j] / (N_averaged + 1.0); - } - - free(D_transform); - free(C_transform); - free(A_transform); - free(B_transform); - free(b_transform); - free(d_transform); - free(c_transform); - free(s_transform); - N_averaged++; - } - - net_free(tmp_net, &c); - - if (save_crit_stress) { - crit_stress[i] = data->extern_field[max_pos]; - } - - if (save_conductivity) { - if (max_pos > 0) { - conductivity[i] = data->conductivity[max_pos - 1]; - } else { - conductivity[i] = cond0; - } - } - - if (save_damage) { - uint_t would_break = 0; - double *tmp_voltage = net_voltages(net, &c); - double *tmp_current = net_currents(net, tmp_voltage, &c); - free(tmp_voltage); - for (uint_t j = 0; j < g->ne; j++) { - bool broken = net->fuses[j]; - bool under_thres = - net->thres[j] < net->thres[data->break_list[max_pos]]; - bool zero_field = fabs(tmp_current[j]) < cutoff; - if (!broken && under_thres && zero_field) { - break_edge(net, j, &c); - } - } - damage[net->num_broken]++; - free(tmp_current); - } - - if (save_energy) { - long double tmp_energy = 0; - if (max_pos > 0) { - long double sigma1 = data->extern_field[0]; - double cond1 = cond0; - for (uint_t j = 0; j < max_pos - 1; j++) { - long double sigma2 = data->extern_field[j + 1]; - double cond2 = data->conductivity[j]; - if (sigma2 > sigma1) { - tmp_energy += 0.5 * gsl_pow_2(sigma1) * (1 - cond2 / cond1) / cond1; - sigma1 = sigma2; - cond1 = cond2; - } - } - } - energy[i] = tmp_energy; - } - - if (save_threshold) { - thresholds[i] = net->thres[data->break_list[max_pos]]; - } - - if (save_current_load) { - loads[i] = - data->extern_field[max_pos] / net->thres[data->break_list[max_pos]]; - } - - if (save_data) { - for (uint_t j = 0; j < data->num_broken; j++) { - fprintf(data_out, "%u %Lg %g ", data->break_list[j], - data->extern_field[j], data->conductivity[j]); - } - fprintf(data_out, "\n"); - } - - data_free(data); - if (save_network) { - FILE *net_out = fopen("network.txt", "w"); - for (uint_t j = 0; j < g->nv; j++) { - fprintf(net_out, "%f %f ", g->vx[2 * j], g->vx[2 * j + 1]); - } - fprintf(net_out, "\n"); - for (uint_t j = 0; j < g->ne; j++) { - fprintf(net_out, "%u %u ", g->ev[2 * j], g->ev[2 * j + 1]); - } - fprintf(net_out, "\n"); - for (uint_t j = 0; j < g->dnv; j++) { - fprintf(net_out, "%f %f ", g->dvx[2 * j], g->dvx[2 * j + 1]); - } - fprintf(net_out, "\n"); - for (uint_t j = 0; j < g->ne; j++) { - fprintf(net_out, "%u %u ", g->dev[2 * j], g->dev[2 * j + 1]); - } - fprintf(net_out, "\n"); - for (uint_t j = 0; j < g->ne; j++) { - fprintf(net_out, "%d ", net->fuses[j]); - } - fclose(net_out); - } - - if (save_cluster_dist) { - uint_t *tmp_cluster_dist = get_cluster_dist(net); - for (uint_t j = 0; j < g->dnv; j++) { - cluster_size_dist[j] += tmp_cluster_dist[j]; - } - free(tmp_cluster_dist); - } - - net_free(net, &c); - graph_free(g); - } - - printf("\033[F\033[JFRACTURE: COMPLETE\n"); - - if (save_cluster_dist) { - FILE *cluster_out = fopen(c_filename, "wb"); - FILE *avalanche_out = fopen(a_filename, "wb"); - - fwrite(cluster_size_dist, sizeof(uint32_t), max_verts, cluster_out); - fwrite(avalanche_size_dist, sizeof(uint32_t), max_edges, avalanche_out); - - fclose(cluster_out); - fclose(avalanche_out); - - free(c_filename); - free(a_filename); - free(cluster_size_dist); - free(avalanche_size_dist); - } - - if (save_conductivity) { - char *cond_filename = (char *)malloc(filename_len * sizeof(char)); - snprintf(cond_filename, filename_len, "cond_%c_%c_%c_%c_%d_%g_%g.dat", - lattice_c, dual_c, boundc, boundc2, L, beta, crack_len); - FILE *cond_file = fopen(cond_filename, "ab"); - fwrite(conductivity, sizeof(double), N, cond_file); - fclose(cond_file); - free(cond_filename); - free(conductivity); - } - - if (save_energy) { - char *tough_filename = (char *)malloc(filename_len * sizeof(char)); - snprintf(tough_filename, filename_len, "enrg_%c_%c_%c_%c_%d_%g_%g.dat", - lattice_c, dual_c, boundc, boundc2, L, beta, crack_len); - FILE *tough_file = fopen(tough_filename, "ab"); - fwrite(energy, sizeof(long double), N, tough_file); - fclose(tough_file); - free(tough_filename); - free(energy); - } - - if (save_threshold) { - char *thres_filename = (char *)malloc(filename_len * sizeof(char)); - snprintf(thres_filename, filename_len, "thrs_%c_%c_%c_%c_%d_%g_%g.dat", - lattice_c, dual_c, boundc, boundc2, L, beta, crack_len); - FILE *thres_file = fopen(thres_filename, "ab"); - fwrite(thresholds, sizeof(long double), N, thres_file); - fclose(thres_file); - free(thres_filename); - free(thresholds); - } - - if (save_current_load) { - char *load_filename = (char *)malloc(filename_len * sizeof(char)); - snprintf(load_filename, filename_len, "load_%c_%c_%c_%c_%d_%g_%g.dat", - lattice_c, dual_c, boundc, boundc2, L, beta, crack_len); - FILE *load_file = fopen(load_filename, "ab"); - fwrite(loads, sizeof(long double), N, load_file); - fclose(load_file); - free(load_filename); - free(loads); - } - - if (save_damage) { - FILE *hdam_file = fopen(d_filename, "wb"); - fwrite(damage, sizeof(uint32_t), max_edges, hdam_file); - fclose(hdam_file); - free(d_filename); - free(damage); - } - - if (save_data) { - fclose(data_out); - } - - if (save_crit_stress) { - char *str_filename = (char *)malloc(filename_len * sizeof(char)); - snprintf(str_filename, filename_len, "strs_%c_%c_%c_%c_%d_%g_%g.dat", - lattice_c, dual_c, boundc, boundc2, L, beta, crack_len); - FILE *str_file = fopen(str_filename, "ab"); - fwrite(crit_stress, sizeof(long double), N, str_file); - fclose(str_file); - free(str_filename); - free(crit_stress); - } - - if (save_correlations) { - FILE *correlations_out = fopen(correlations_filename, "wb"); - fwrite(&N_averaged, sizeof(uint64_t), 1, correlations_out); - fwrite(&mean_d, sizeof(double), 1, correlations_out); - fwrite(&mean_c, sizeof(double), 1, correlations_out); - fwrite(&mean_b, sizeof(double), 1, correlations_out); - fwrite(&mean_s, sizeof(double), 1, correlations_out); - fwrite(&mean_A, sizeof(double), 1, correlations_out); - fwrite(&mean_D, sizeof(double), 1, correlations_out); - fwrite(&mean_C, sizeof(double), 1, correlations_out); - fwrite(&mean_B, sizeof(double), 1, correlations_out); - fwrite(dd_correlations, sizeof(double), pow(L, 2), correlations_out); - fwrite(dc_correlations, sizeof(double), pow(L, 2), correlations_out); - fwrite(db_correlations, sizeof(double), pow(L, 2), correlations_out); - fwrite(ds_correlations, sizeof(double), pow(L, 2), correlations_out); - fwrite(dA_correlations, sizeof(double), pow(L, 2), correlations_out); - fwrite(cc_correlations, sizeof(double), pow(L, 2), correlations_out); - fwrite(cb_correlations, sizeof(double), pow(L, 2), correlations_out); - fwrite(cs_correlations, sizeof(double), pow(L, 2), correlations_out); - fwrite(cA_correlations, sizeof(double), pow(L, 2), correlations_out); - fwrite(bb_correlations, sizeof(double), pow(L, 2), correlations_out); - fwrite(bs_correlations, sizeof(double), pow(L, 2), correlations_out); - fwrite(bA_correlations, sizeof(double), pow(L, 2), correlations_out); - fwrite(ss_correlations, sizeof(double), pow(L, 2), correlations_out); - fwrite(AA_correlations, sizeof(double), pow(L, 2), correlations_out); - fwrite(DD_correlations, sizeof(double), pow(L, 2), correlations_out); - fwrite(DC_correlations, sizeof(double), pow(L, 2), correlations_out); - fwrite(DB_correlations, sizeof(double), pow(L, 2), correlations_out); - fwrite(CC_correlations, sizeof(double), pow(L, 2), correlations_out); - fwrite(CB_correlations, sizeof(double), pow(L, 2), correlations_out); - fwrite(BB_correlations, sizeof(double), pow(L, 2), correlations_out); - fclose(correlations_out); - - free(dd_correlations); - free(dc_correlations); - free(db_correlations); - free(ds_correlations); - free(dA_correlations); - free(cc_correlations); - free(cb_correlations); - free(cs_correlations); - free(cA_correlations); - free(bb_correlations); - free(bs_correlations); - free(bA_correlations); - free(ss_correlations); - free(AA_correlations); - free(DD_correlations); - free(DC_correlations); - free(DB_correlations); - free(CC_correlations); - free(CB_correlations); - free(BB_correlations); - fftw_free(fftw_forward_in); - fftw_free(fftw_forward_out); - fftw_free(fftw_reverse_in); - fftw_free(fftw_reverse_out); - fftw_destroy_plan(forward_plan); - fftw_destroy_plan(reverse_plan); - free(correlations_filename); - } - - fftw_cleanup(); - - CHOL_F(finish)(&c); - - return 0; -} diff --git a/src/fracture.cpp b/src/fracture.cpp new file mode 100644 index 0000000..77af253 --- /dev/null +++ b/src/fracture.cpp @@ -0,0 +1,59 @@ + +#include <random> +#include <iostream> + +#include <cholmod.h> + +#include "randutils/randutils.hpp" + +#include <graph.hpp> +#include <network.hpp> +#include <hooks.hpp> +#include "measurements.hpp" + +int main(int argc, char* argv[]) { + int opt; + + unsigned int N = 1; + unsigned int L = 16; + double beta = 0.5; + + while ((opt = getopt(argc, argv, "N:L:b:")) != -1) { + switch (opt) { + case 'N': + N = (unsigned int)atof(optarg); + break; + case 'L': + L = atoi(optarg); + break; + case 'b': + beta = atof(optarg); + break; + default: + exit(1); + } + } + + cholmod_common c; + CHOL_F(start)(&c); + c.supernodal = CHOLMOD_SUPERNODAL; + + + graph G(L); + network base_network(G, &c); + ma meas(N, L, beta); + + randutils::auto_seed_128 seeds; + std::mt19937 rng{seeds}; + + for (unsigned int trial = 0; trial < N; trial++) { + network tmp_network(base_network); + tmp_network.set_thresholds(beta, rng); + tmp_network.fracture(meas); + } + + CHOL_F(finish)(&c); + + return 0; +} + diff --git a/src/fracture.h b/src/fracture.h deleted file mode 100644 index 5eb0a1d..0000000 --- a/src/fracture.h +++ /dev/null @@ -1,123 +0,0 @@ - -#pragma once - -#include <assert.h> -#include <cholmod.h> -#include <float.h> -#include <getopt.h> -#include <gsl/gsl_math.h> -#include <gsl/gsl_randist.h> -#include <gsl/gsl_rng.h> -#include <gsl/gsl_sf_erf.h> -#include <gsl/gsl_sf_exp.h> -#include <gsl/gsl_sf_log.h> -#include <inttypes.h> -#include <math.h> -#include <stdbool.h> -#include <stdint.h> -#include <stdio.h> -#include <stdlib.h> -#include <string.h> -#include <sys/types.h> -#include <unistd.h> - -#include <jst/graph.h> -#include <jst/rand.h> - -// these defs allow me to switch to long int cholmod in a sitch -#define int_t int -#define uint_t unsigned int -#define CINT_MAX INT_MAX -#define CHOL_F(x) cholmod_##x - -#define GSL_RAND_GEN gsl_rng_mt19937 - -typedef struct { - const graph_t *graph; - bool *fuses; - long double *thres; - double inf; - cholmod_dense *boundary_cond; - cholmod_factor *factor; - bool voltage_bound; - uint_t num_broken; - uint_t dim; - uint_t nep; - uint_t *evp; - cholmod_sparse *voltcurmat; -} net_t; - -typedef struct { - uint_t num_broken; - uint_t *break_list; - double *conductivity; - long double *extern_field; -} data_t; - -intptr_t *run_voronoi(uint_t num_coords, double *coords, bool periodic, - double xmin, double xmax, double ymin, double ymax); - -cholmod_sparse *gen_adjacency(const net_t *net, bool dual, bool use_gp, - bool symmetric, cholmod_common *c); - -cholmod_sparse *gen_laplacian(const net_t *net, cholmod_common *c); - -int edge_to_verts(uint_t width, bool periodic, uint_t edge, bool index); - -int dual_edge_to_verts(uint_t width, bool periodic, uint_t edge, bool index); - -double dual_vert_to_coord(uint_t width, bool periodic, uint_t vert, bool index); - -void factor_update(cholmod_factor *factor, uint_t v1, uint_t v2, - cholmod_common *c); -void factor_update2(cholmod_factor *factor, uint_t v1, cholmod_common *c); - -void net_notch(net_t *net, double notch_len, cholmod_common *c); -data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff); -double *net_voltages(const net_t *net, cholmod_common *c); -double *net_currents(const net_t *net, const double *voltages, - cholmod_common *c); -double net_conductivity(const net_t *net, const double *voltages); - -void update_boundary(net_t *instance, const double *avg_field); - -FILE *get_file(const char *prefix, uint_t width, uint_t crack, double beta, - uint_t iter, uint_t num_iter, uint_t num, bool read); - -double update_beta(double beta, uint_t width, const double *stress, - const double *damage, double bound_total); - -cholmod_sparse *gen_voltcurmat(uint_t num_edges, uint_t num_verts, - uint_t *edges_to_verts, cholmod_common *c); - -net_t *net_copy(const net_t *net, cholmod_common *c); - -void net_free(net_t *instance, cholmod_common *c); - -net_t *net_create(const graph_t *g, double inf, double beta, double notch_len, - bool vb, cholmod_common *c); - -bool break_edge(net_t *instance, uint_t edge, cholmod_common *c); - -components_t *get_clusters(net_t *instance); - -uint_t *get_cluster_dist(net_t *instance); - -void randfunc_flat(gsl_rng *r, double *x, double *y); -void randfunc_gaus(gsl_rng *r, double *x, double *y); - -double *get_corr(net_t *instance, uint_t **dists, cholmod_common *c); - -double *bin_values(graph_t *network, uint_t width, double *values); - -cholmod_dense *bound_set(const graph_t *g, bool vb, double notch_len, - cholmod_common *c); - -data_t *data_create(uint_t num_edges); -void data_free(data_t *data); -void data_update(data_t *data, uint_t last_broke, long double strength, - double conductivity); - -long double rand_dist_pow(const gsl_rng *r, double beta); - -bool is_in(uint_t len, uint_t *list, uint_t element); diff --git a/src/long_anal_process.c b/src/long_anal_process.c deleted file mode 100644 index d4ec4e0..0000000 --- a/src/long_anal_process.c +++ /dev/null @@ -1,158 +0,0 @@ - -#include "fracture.h" -#include <gsl/gsl_blas.h> -#include <gsl/gsl_matrix.h> -#include <gsl/gsl_sf_erf.h> -#include <gsl/gsl_sf_laguerre.h> -#include <gsl/gsl_vector.h> -#include <sys/stat.h> - -void get_mean(uint_t len, long double *data, long double result[2]) { - long double total = 0; - - for (uint_t i = 0; i < len; i++) { - total += data[i]; - } - - long double mean = total / len; - long double total_sq_diff = 0; - - for (uint_t i = 0; i < len; i++) { - total_sq_diff += pow(mean - data[i], 2); - } - - long double mean_err = sqrt(total_sq_diff) / len; - - result[0] = mean; - result[1] = mean_err; -} - -void get_mom(uint_t len, uint_t n, long double *data, long double mean[2], - long double result[2]) { - long double total_n_diff = 0; - long double total_n2_diff = 0; - - for (uint_t i = 0; i < len; i++) { - total_n_diff += pow(fabsl(mean[0] - data[i]), n); - total_n2_diff += pow(fabsl(mean[0] - data[i]), 2 * n); - } - - long double mom = total_n_diff / len; - long double mom_err = sqrt(total_n2_diff) / len; - - result[0] = mom; - result[1] = mom_err; -} - -int main(int argc, char *argv[]) { - uint_t nc = argc - 1; - uint_t *Ls_c = (uint_t *)malloc(nc * sizeof(uint_t)); - double *betas_c = (double *)malloc(nc * sizeof(double)); - long double *vals_c1 = (long double *)malloc(nc * sizeof(long double)); - long double *errors_c1 = (long double *)malloc(nc * sizeof(long double)); - long double *vals_c2 = (long double *)malloc(nc * sizeof(long double)); - long double *errors_c2 = (long double *)malloc(nc * sizeof(long double)); - long double *vals_c3 = (long double *)malloc(nc * sizeof(long double)); - long double *errors_c3 = (long double *)malloc(nc * sizeof(long double)); - long double *vals_c4 = (long double *)malloc(nc * sizeof(long double)); - long double *errors_c4 = (long double *)malloc(nc * sizeof(long double)); - - char *out_filename = (char *)malloc(100 * sizeof(char)); - uint_t out_filename_len = 0; - - for (uint_t i = 0; i < nc; i++) { - char *fn = argv[1 + i]; - char *buff = (char *)malloc(20 * sizeof(char)); - uint_t pos = 0; - uint_t c = 0; - while (c < 5) { - if (fn[pos] == '_') { - c++; - } - if (i == 0) { - out_filename[pos] = fn[pos]; - out_filename_len++; - } - pos++; - } - c = 0; - while (fn[pos] != '_') { - buff[c] = fn[pos]; - pos++; - c++; - } - buff[c] = '\0'; - Ls_c[i] = atoi(buff); - c = 0; - pos++; - while (fn[pos] != '_') { - buff[c] = fn[pos]; - pos++; - c++; - } - buff[c] = '\0'; - betas_c[i] = atof(buff); - free(buff); - - struct stat info; - stat(fn, &info); - uint_t num_bytes = info.st_size; - uint_t dist_len = (num_bytes * sizeof(char)) / sizeof(long double); - - long double *dist = malloc(dist_len * sizeof(long double)); - FILE *dist_file = fopen(fn, "rb"); - fread(dist, sizeof(long double), dist_len, dist_file); - fclose(dist_file); - { - long double mom1[2]; - get_mean(dist_len, dist, mom1); - vals_c1[i] = mom1[0]; - errors_c1[i] = mom1[1]; - - long double mom2[2]; - get_mom(dist_len, 2, dist, mom1, mom2); - vals_c2[i] = mom2[0]; - errors_c2[i] = mom2[1]; - - long double mom3[2]; - get_mom(dist_len, 3, dist, mom1, mom3); - vals_c3[i] = mom3[0]; - errors_c3[i] = mom3[1]; - - long double mom4[2]; - get_mom(dist_len, 4, dist, mom1, mom4); - vals_c4[i] = mom4[0]; - errors_c4[i] = mom4[1]; - } - free(dist); - } - - out_filename[out_filename_len - 1] = '.'; - out_filename[out_filename_len] = 't'; - out_filename[out_filename_len + 1] = 'x'; - out_filename[out_filename_len + 2] = 't'; - out_filename[out_filename_len + 3] = '\0'; - - FILE *out_file = fopen(out_filename, "w"); - - for (uint_t i = 0; i < nc; i++) { - fprintf(out_file, "%u %g %Lg %Lg %Lg %Lg %Lg %Lg %Lg %Lg\n", Ls_c[i], - betas_c[i], vals_c1[i], errors_c1[i], vals_c2[i], errors_c2[i], - vals_c3[i], errors_c3[i], vals_c4[i], errors_c4[i]); - } - - fclose(out_file); - - free(Ls_c); - free(betas_c); - free(vals_c1); - free(errors_c1); - free(vals_c2); - free(errors_c2); - free(vals_c3); - free(errors_c3); - free(vals_c4); - free(errors_c4); - - return 0; -} diff --git a/src/measurements.cpp b/src/measurements.cpp new file mode 100644 index 0000000..7d5d539 --- /dev/null +++ b/src/measurements.cpp @@ -0,0 +1,248 @@ + +#include "measurements.hpp" + +bool trivial(boost::detail::edge_desc_impl<boost::undirected_tag,unsigned long>) { + return true; +} + +void update_distribution_data(std::string id, const std::vector<unsigned int>& data, unsigned int N, unsigned int L, double beta) { + std::string filename = "fracture_" + std::to_string(L) + "_" + std::to_string(beta) + "_" + id + ".dat"; + std::ifstream file(filename); + + uint64_t N_old = 0; + std::vector<uint64_t> data_old(data.size(), 0); + + if (file.is_open()) { + file >> N_old; + for (unsigned int i = 0; i < data.size(); i++) { + uint64_t num; + file >> num; + data_old[i] = num; + } + + file.close(); + } + + std::ofstream file_out(filename); + + file_out <<std::fixed<< N_old + N << "\n"; + for (unsigned int i = 0; i < data.size(); i++) { + file_out <<std::fixed<< data_old[i] + data[i] << " "; + } + + file_out.close(); +} + +void update_field_data(std::string id, double tot, const std::vector<double>& data, unsigned int L, double beta) { + std::string filename = "fracture_" + std::to_string(L) + "_" + std::to_string(beta) + "_" + id + ".dat"; + std::ifstream file(filename); + + uint64_t N_old = 0; + uint64_t tot_old = 0; + uint64_t tot_2_old = 0; + std::vector<uint64_t> data_old(data.size(), 0); + std::vector<uint64_t> data_old_2(data.size(), 0); + + if (file.is_open()) { + file >> N_old; + file >> tot_old; + file >> tot_2_old; + for (unsigned int i = 0; i < data.size(); i++) { + uint64_t num; + file >> num; + data_old[i] = num; + } + for (unsigned int i = 0; i < data.size(); i++) { + uint64_t num; + file >> num; + data_old_2[i] = num; + } + + file.close(); + } + + std::ofstream file_out(filename); + + file_out <<std::fixed<< N_old + 1 << " " << (uint64_t)(tot_old + tot) << " " << (uint64_t)(tot_2_old + pow(tot, 2))<< "\n"; + for (unsigned int i = 0; i < data.size(); i++) { + file_out << (uint64_t)(data_old[i] + data[i]) << " "; + } + file_out <<std::fixed<< "\n"; + for (unsigned int i = 0; i < data.size(); i++) { + file_out << (uint64_t)(data_old_2[i] + pow(data[i], 2)) << " "; + } + + file_out.close(); +} + +std::vector<fftw_complex> data_transform(unsigned int L, const std::vector<double>& data, fftw_plan forward_plan, double *fftw_forward_in, fftw_complex *fftw_forward_out) { + for (unsigned int i = 0; i < pow(L, 2); i++) { + fftw_forward_in[i] = data[i]; + } + + fftw_execute(forward_plan); + + std::vector<fftw_complex> output(pow(L, 2)); + + for (unsigned int i = 0; i < pow(L, 2); i++) { + output[i][0] = fftw_forward_out[i][0]; + output[i][1] = fftw_forward_out[i][1]; + } + + return output; +} + +std::vector<double> correlation(unsigned int L, const std::vector<fftw_complex>& tx1, const std::vector<fftw_complex>& tx2, fftw_plan reverse_plan, fftw_complex *fftw_reverse_in, double *fftw_reverse_out) { + for (unsigned int i = 0; i < pow(L, 2); i++) { + fftw_reverse_in[i][0] = tx1[i][0] * tx2[i][0] + tx1[i][1] * tx2[i][1]; + fftw_reverse_in[i][1] = tx1[i][0] * tx2[i][1] - tx1[i][1] * tx2[i][0]; + } + + fftw_execute(reverse_plan); + + std::vector<double> output(pow(L / 2 + 1, 2)); + + for (unsigned int i = 0; i < pow(L / 2 + 1, 2); i++) { + output[i] = fftw_reverse_out[L * (i / (L / 2 + 1)) + i % (L / 2 + 1)] / pow(L, 2); + } + + return output; +} + +ma::ma(unsigned int N, unsigned int L, double beta) : L(L), G(2 * pow(L / 2, 2)), bin_counts(log2(L) + 1, 0), N(N), beta(beta) { + ad.resize(pow(L, 2), 0); + cd.resize(pow(L, 2), 0); + + // FFTW setup for correlation functions + fftw_set_timelimit(1); + + fftw_forward_in = (double *)fftw_malloc(pow(L, 2) * sizeof(double)); + fftw_forward_out = (fftw_complex *)fftw_malloc(pow(L, 2) * sizeof(fftw_complex)); + fftw_reverse_in = (fftw_complex *)fftw_malloc(pow(L, 2) * sizeof(fftw_complex)); + fftw_reverse_out = (double *)fftw_malloc(pow(L, 2) * sizeof(double)); + + forward_plan = fftw_plan_dft_r2c_2d(L, L, fftw_forward_in, fftw_forward_out, 0); + reverse_plan = fftw_plan_dft_c2r_2d(L, L, fftw_reverse_in, fftw_reverse_out, 0); +} + +ma::~ma() { + // clean up FFTW objects + fftw_free(fftw_forward_in); + fftw_free(fftw_forward_out); + fftw_free(fftw_reverse_in); + fftw_free(fftw_reverse_out); + fftw_destroy_plan(forward_plan); + fftw_destroy_plan(reverse_plan); + fftw_cleanup(); + + update_distribution_data("na", ad, N, L, beta); + update_distribution_data("nc", cd, N, L, beta); + update_distribution_data("bc", bin_counts, N, L, beta); + +} + +void ma::pre_fracture(const network &) { + lv = 0; + as = 0; + avalanches.clear(); + boost::remove_edge_if(trivial, G); +} + +void ma::bond_broken(const network& net, const std::pair<double, std::vector<double>>& cur, unsigned int i) { + if (cur.first / fabs(cur.second[i]) * net.thresholds[i] > lv) { + ad[as]++; + as = 0; + lv = cur.first / fabs(cur.second[i]) * net.thresholds[i]; + avalanches.push_back({}); + } else { + as++; + avalanches.back().push_back(i); + } + + boost::add_edge(net.G.dual_edges[i][0], net.G.dual_edges[i][1], G); +} + +void ma::post_fracture(network &n) { + std::vector<unsigned int> component(boost::num_vertices(G)); + unsigned int num = boost::connected_components(G, &component[0]); + + std::vector<unsigned int> comp_sizes(num, 0); + + for (unsigned int c : component) { + comp_sizes[c]++; + } + + unsigned int max_i = 0; + unsigned int max_size = 0; + + for (unsigned int i = 0; i < num; i++) { + if (comp_sizes[i] > max_size) { + max_i = i; + max_size = comp_sizes[i]; + } + } + + for (unsigned int be = 0; be < log2(L); be++) { + unsigned int bin = pow(2, be); + + for (unsigned int i = 0; i < pow(L / bin, 2); i++) { + bool in_bin = false; + for (unsigned int j = 0; j < pow(bin, 2); j++) { + unsigned int edge = L * (bin * ((i * bin) / L) + j / bin) + (i * bin) % L + j % bin; + if (!n.fuses[edge] && max_i == component[n.G.dual_edges[edge][0]]) { + in_bin = true; + break; + } + } + + if (in_bin) { + bin_counts[be]++; + } + } + } + + bin_counts[log2(L)]++; + + std::vector<double> crack_damage(pow(L, 2), 0.0); + + double damage_tot = 0; + for (unsigned int i = 0; i < pow(L, 2); i++) { + if (!n.fuses[i] && max_i == component[n.G.dual_edges[i][0]]) { + damage_tot++; + crack_damage[i] = 1.0; + } + } + + std::vector<fftw_complex> t_crack_damage = data_transform(L, crack_damage, forward_plan, fftw_forward_in, fftw_forward_out); + std::vector<double> Ccc = correlation(L, t_crack_damage, t_crack_damage, reverse_plan, fftw_reverse_in, fftw_reverse_out); + + update_field_data("Ccc", damage_tot, Ccc, L, beta); + + for (auto e = avalanches.back().rbegin(); e != avalanches.back().rend(); e++) { + boost::remove_edge(n.G.dual_edges[*e][0], n.G.dual_edges[*e][1], G); + } + + num = boost::connected_components(G, &component[0]); + + for (unsigned int i = 0; i < num; i++) { + double size = 0; + std::fill(crack_damage.begin(), crack_damage.end(), 0.0); + + for (unsigned int j = 0; j < pow(L, 2); j++) { + if (component[n.G.edges[j][0]] == i && n.fuses[j]) { + size++; + crack_damage[j] = 1.0; + } + } + + if (size > 0) { + cd[size - 1]++; + t_crack_damage = data_transform(L, crack_damage, forward_plan, fftw_forward_in, fftw_forward_out); + Ccc = correlation(L, t_crack_damage, t_crack_damage, reverse_plan, fftw_reverse_in, fftw_reverse_out); + + update_field_data("Cclcl", size, Ccc, L, beta); + } + } + +} + diff --git a/src/measurements.hpp b/src/measurements.hpp new file mode 100644 index 0000000..817481b --- /dev/null +++ b/src/measurements.hpp @@ -0,0 +1,54 @@ + +#include <vector> +#include <algorithm> +#include <cmath> +#include <list> +#include <fstream> +#include <string> +#include <cinttypes> +#include <sstream> + +#include <boost/graph/adjacency_list.hpp> +#include <boost/graph/connected_components.hpp> + +#include <fftw3.h> + +#include <network.hpp> +#include <hooks.hpp> + +class ma : public hooks { + // need: + // - measurements for correlation functions (nice looking?) + // - interface for reading and writing from files (nice looking?) + // - interface for turning on and off specific measurements + // + private: + double *fftw_forward_in; + fftw_complex *fftw_forward_out; + fftw_complex *fftw_reverse_in; + double *fftw_reverse_out; + fftw_plan forward_plan; + fftw_plan reverse_plan; + boost::adjacency_list <boost::listS, boost::vecS, boost::undirectedS> G; + + public: + unsigned int N; + unsigned int L; + unsigned int as; + double lv; + double beta; + + + std::vector<unsigned int> ad; + std::vector<unsigned int> cd; + std::vector<unsigned int> bin_counts; + std::list<std::list<unsigned int>> avalanches;; + + ma(unsigned int N, unsigned int L, double beta); + ~ma(); + + void pre_fracture(const network &) override; + void bond_broken(const network& net, const std::pair<double, std::vector<double>>& cur, unsigned int i) override; + void post_fracture(network &n) override; +}; + diff --git a/src/randutils b/src/randutils new file mode 160000 +Subproject 8486a610a954a8248c12485fb4cfc390a5f5f85 |