Dear Lorenzo, Thanks. Regarding Python and C++, there are not so much GetFEM specific differences between the two. Once you have a model implemented and working in Python it should be minor work to port it to C++ for performance reasons if you are familiar with the language. I agree that we should emphasize the use of GWFL in the documentation but I also highly recommend to read our paper on GWFL <https://www.researchgate.net/publication/347772050_GetFEM_Automated_FE_Modeling_of_Multiphysics_Problems_Based_on_a_Generic_Weak_Form_Language> before trying to use GetFEM.
A. the error message that you get is consistent with the documentation [image: image.png] the value of 0 is excluded. Actually, as far as I can see, for theta=0 you do not have an implicit method any more. Personally I prefer to implement time integration manually than using convenience functions such as "add_theta_method_for_first_order". Such functions add variables and data to your model for you automatically, but this black-box mode is probably not the best for learning. B. If the model only has linear terms, the Newton solver is skipped, and model.solve() will just perform a single linear solve with whatever linear solver you have specified. C. If you only do model.assembly(option='build_rhs') the tangent matrix is not assembled, so in a subsequent solve it will be necessary to assemble anyway. I do not remember the details, but the general idea is that you should not mix manual model.assembly() and model.solve() calls. If you use model.solve() it will do only the minimum amount of assemblies necessary, based on detecting changes in variables and data since the last assembly. You only use model.assembly() in combination with e.g. "gf.linsolve_mumps(matrix, residual)" if you implement your own Netwon solver in Python, in that case you skip md.solve(). In this case you take care yourself about performing only the necessary assemblies using the options ‘build_all’, ‘build_rhs’, ‘build_matrix’. According to my experience for a nonlinear problem it doesn't matter so much if you have only one nonlinear or several nonlinear terms. If you have to do 6-7 linear solves per simulation step anyway, it doesn't matter so much how many terms you have to reassemble. If you find, for some specific problem, that assembly time is your bottle-neck compared to the time for linear solves, we can have a look at it and help you to optimize it. Best regards Kostas On Fri, Oct 22, 2021 at 6:27 PM Lorenzo Ferro <lorenzo.i...@gmail.com> wrote: > Dear Poulios and Yves, > > thank you very much for your kind and prompt replies. > > Little by little I'm getting familiar with the logics of this FEM library > that I found very interesting and I'd like to go in depth more in the near > future, learning also the C++ features, but it will take a good amount of > time that in this moment I don't have, do for now I will continue mastering > the Python interface. > > I believe that you should emphasize in the manuals and push the usage of > the GWFL in place of all the other model bricks, since for the beginners it > is not obvious and indeed it makes life easier using "add_linear_term" and > "add_nonlinear_term" for everything. Thank you for this enlightening > clarification. > > I've some additional questions that rose up today: > > A) If I initialize a theta-method with theta=0 (explicit forward euler) I > get the following error > RuntimeError: (Getfem::InterfaceError) -- Error in getfem_models.cc, line > 1281 : > > Invalid value of theta parameter for the theta-method > > but the manual doesn't exclude this option (even though the manual is not > complete at the bottom of page The model tools for the integration of > transient problems — GetFEM > <https://getfem.org/userdoc/model_time_integration.html>) > > > B) Let's imagine that our model is perfectly linear and I assemble it with > "model.assembly()". What happens if I run "model.solve()" ? Are the > equations assembled again or the software somehow understand that the > assembly has already been done and skips to the proper solver? > > C) Linked to the previous question. If I need to change only the RHS by > "model.assembly(option='build_rhs')" and after that I run "model.solve()", > does this method assembly again the whole model including the tangent > matrix or on the contrary it is smart and skips to the solver? > > I ask these things because I'm trying to generate a 3D thermal-electric > fully coupled model that I found quite slow in solving the equations since > the model is strongly non-linear everywhere. To date, it is still not fully > featured. > In the near future I'll probably make some sensitivity analysis so as to > try to linearize the model where possible, and for this reason I'm trying > to understand possible performance improvements using the built-in standard > solver. > > Thank you in advance for your kind feedback. > ...and compliments for the nice FEM library you're developing. > > HAve a nice week-end > Lorenzo > > > Il giorno ven 22 ott 2021 alle ore 10:43 Konstantinos Poulios < > logar...@googlemail.com> ha scritto: > >> Dear Lorenzo, >> >> Thanks for your interest. Below some answers to your questions >> >> BR >> Kostas >> >> On Thu, Oct 21, 2021 at 12:47 AM Lorenzo Ferro <lorenzo.i...@gmail.com> >> wrote: >> >>> Dear All, >>> >>> I've started testing the GeFEM library and I've some questions, with >>> particular reference to the Python interface (the only one I can use so >>> far): >>> >>> 1) the "add_lumped_mass_for_first_order_brick" method is not available >>> even though it is mentioned in the Python-Model guide for GetFEM 5.4. I've >>> installed the precompiled 5.4 version for Windows with >>> getfem5.4win-amd64-py3.7.exe. How can I get the Lumped_Mass feature for the >>> Python interface? >>> >> The scripting language interface for this function has been added after >> the last release: >> >> https://git.savannah.nongnu.org/cgit/getfem.git/commit/interface/src/gf_model_set.cc?id=0241496a72624aca5f7f6a4f3725accee863d46e >> so you will have to wait for the next release or compile GetFEM yourself. >> >> >>> 2) The documentation of bricks like "Generic Elliptic Brick", >>> "Fourier_Robin Brick" and "Source Term Brick" states that the coefficients >>> could be arbitrary GWFL expressions. What happens if I'd input non-linear >>> coefficients? How the software manages the assembly in this case? Or, it >>> would be better to use the Model method "add_non_linear_term"? >>> >> Just use the add_linear_term or add_nonlinear_term for this. E.g. The >> add_fourier_robin_brick is just a convenience function calling >> add_(non)linear_term for you >> [image: image.png] >> Such bricks made sense before the introduction of GWFL, nowadays you >> should just use add_linear_term or add_nonlinear_term for everything. >> >> 3) If I'd like to implement a radiation non-linear boundary condition, >>> what would be the more advisable way? Is in this case the >>> "add_non_linear_term" method the correct way? >>> >> Yes. If you have a specific equation you want to implement, we can help >> with writing it in GWFL. >> >> >>> 4) I've not understood what is the difference between "add_source_term" >>> and "add_source_term_brick". Can you advise when I should use one or the >>> other (with a practical example). >>> >> If I see correctly, providing an expression "expr" and a variable name >> "var" to the second one is equivalent to providing "expr*Test_var" in the >> first one. In general you should not use these bricks because they add a >> contribution to the rhs that is not accounted for in the stiffness matrix >> computation. Rarely do you need something like this. You can just use >> "add_linear_term" and "add_nonlinear_term" for everything. >> >> >>> 5) Is there a way to visualize mesh regions for doublecheck purposes? >>> Maybe, is there a way to export such regions by means of the Python >>> interface? >>> >> Good question. There are different ways, maybe the most simple is >> something like >> mf.export_to_vtk("meshregions.vtk", mf, md.interpolation("1",mf,1), >> 'region1', mf, md.interpolation("1",mf,2), 'region2', mf, >> md.interpolation("1",mf,4), 'region4') >> where mf is a (preferably discontinuous) mesh_fem and md is just a model >> (it can be an empty one). >> >> >>> Thank you in advance for your kind support. >>> Regards, >>> Lorenzo >>> >>
