summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorJaron Kent-Dobias <jaron@kent-dobias.com>2018-11-01 12:33:37 -0400
committerJaron Kent-Dobias <jaron@kent-dobias.com>2018-11-01 12:33:37 -0400
commit07906baa42470bad14d2c40f57967691f6118969 (patch)
tree416ae624829967861c7c799103b3ff795e9e36b4
parent8c4c42d81745ea33c31150fe22e834d97e29ede6 (diff)
downloadfuse_networks-07906baa42470bad14d2c40f57967691f6118969.tar.gz
fuse_networks-07906baa42470bad14d2c40f57967691f6118969.tar.bz2
fuse_networks-07906baa42470bad14d2c40f57967691f6118969.zip
revamped and simplied fracture code with c++
-rw-r--r--.gitmodules3
-rw-r--r--CMakeLists.txt22
-rw-r--r--LICENSE674
-rw-r--r--README.md4
-rw-r--r--lib/CMakeLists.txt4
-rw-r--r--lib/bound_set.c110
-rw-r--r--lib/break_edge.c50
-rw-r--r--lib/correlations.c46
-rw-r--r--lib/data.c35
-rw-r--r--lib/factor_update.c69
-rw-r--r--lib/fracture.h123
-rw-r--r--lib/gen_laplacian.c142
-rw-r--r--lib/gen_voltcurmat.c36
-rw-r--r--lib/geometry.c55
-rw-r--r--lib/get_dual_clusters.c31
-rw-r--r--lib/include/graph.hpp29
-rw-r--r--lib/include/hooks.hpp12
-rw-r--r--lib/include/network.hpp48
-rw-r--r--lib/net.c146
-rw-r--r--lib/net_conductivity.c35
-rw-r--r--lib/net_currents.c52
-rw-r--r--lib/net_fracture.c68
-rw-r--r--lib/net_voltages.c39
-rw-r--r--lib/rand.c15
-rw-r--r--lib/src/graph.cpp38
-rw-r--r--lib/src/network.cpp280
-rw-r--r--src/CMakeLists.txt7
-rw-r--r--src/anal_process.c135
-rw-r--r--src/big_anal_process.c158
-rw-r--r--src/cen_anal_process.c157
-rw-r--r--src/corr_test.c27
-rw-r--r--src/fracture.c1068
-rw-r--r--src/fracture.cpp59
-rw-r--r--src/fracture.h123
-rw-r--r--src/long_anal_process.c158
-rw-r--r--src/measurements.cpp248
-rw-r--r--src/measurements.hpp54
m---------src/randutils0
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