You can subscribe to this list here.
2003 
_{Jan}

_{Feb}

_{Mar}

_{Apr}

_{May}

_{Jun}

_{Jul}

_{Aug}

_{Sep}
(2) 
_{Oct}
(2) 
_{Nov}
(27) 
_{Dec}
(31) 

2004 
_{Jan}
(6) 
_{Feb}
(15) 
_{Mar}
(33) 
_{Apr}
(10) 
_{May}
(46) 
_{Jun}
(11) 
_{Jul}
(21) 
_{Aug}
(15) 
_{Sep}
(13) 
_{Oct}
(23) 
_{Nov}
(1) 
_{Dec}
(8) 
2005 
_{Jan}
(27) 
_{Feb}
(57) 
_{Mar}
(86) 
_{Apr}
(23) 
_{May}
(37) 
_{Jun}
(34) 
_{Jul}
(24) 
_{Aug}
(17) 
_{Sep}
(50) 
_{Oct}
(24) 
_{Nov}
(10) 
_{Dec}
(60) 
2006 
_{Jan}
(47) 
_{Feb}
(46) 
_{Mar}
(127) 
_{Apr}
(19) 
_{May}
(26) 
_{Jun}
(62) 
_{Jul}
(47) 
_{Aug}
(51) 
_{Sep}
(61) 
_{Oct}
(42) 
_{Nov}
(50) 
_{Dec}
(33) 
2007 
_{Jan}
(60) 
_{Feb}
(55) 
_{Mar}
(77) 
_{Apr}
(102) 
_{May}
(82) 
_{Jun}
(102) 
_{Jul}
(169) 
_{Aug}
(117) 
_{Sep}
(80) 
_{Oct}
(37) 
_{Nov}
(51) 
_{Dec}
(43) 
2008 
_{Jan}
(71) 
_{Feb}
(94) 
_{Mar}
(98) 
_{Apr}
(125) 
_{May}
(54) 
_{Jun}
(119) 
_{Jul}
(60) 
_{Aug}
(111) 
_{Sep}
(118) 
_{Oct}
(125) 
_{Nov}
(119) 
_{Dec}
(94) 
2009 
_{Jan}
(109) 
_{Feb}
(38) 
_{Mar}
(93) 
_{Apr}
(88) 
_{May}
(29) 
_{Jun}
(57) 
_{Jul}
(53) 
_{Aug}
(48) 
_{Sep}
(68) 
_{Oct}
(151) 
_{Nov}
(23) 
_{Dec}
(35) 
2010 
_{Jan}
(84) 
_{Feb}
(60) 
_{Mar}
(184) 
_{Apr}
(112) 
_{May}
(60) 
_{Jun}
(90) 
_{Jul}
(23) 
_{Aug}
(70) 
_{Sep}
(119) 
_{Oct}
(27) 
_{Nov}
(47) 
_{Dec}
(54) 
2011 
_{Jan}
(22) 
_{Feb}
(19) 
_{Mar}
(92) 
_{Apr}
(93) 
_{May}
(35) 
_{Jun}
(91) 
_{Jul}
(32) 
_{Aug}
(61) 
_{Sep}
(7) 
_{Oct}
(69) 
_{Nov}
(81) 
_{Dec}
(23) 
2012 
_{Jan}
(64) 
_{Feb}
(95) 
_{Mar}
(35) 
_{Apr}
(36) 
_{May}
(63) 
_{Jun}
(98) 
_{Jul}
(70) 
_{Aug}
(171) 
_{Sep}
(149) 
_{Oct}
(64) 
_{Nov}
(67) 
_{Dec}
(126) 
2013 
_{Jan}
(108) 
_{Feb}
(104) 
_{Mar}
(171) 
_{Apr}
(133) 
_{May}
(108) 
_{Jun}
(100) 
_{Jul}
(93) 
_{Aug}
(126) 
_{Sep}
(74) 
_{Oct}
(59) 
_{Nov}
(145) 
_{Dec}
(93) 
2014 
_{Jan}
(38) 
_{Feb}
(45) 
_{Mar}
(26) 
_{Apr}
(41) 
_{May}
(125) 
_{Jun}
(70) 
_{Jul}
(61) 
_{Aug}
(66) 
_{Sep}
(60) 
_{Oct}
(110) 
_{Nov}
(27) 
_{Dec}
(30) 
2015 
_{Jan}
(43) 
_{Feb}
(67) 
_{Mar}
(71) 
_{Apr}
(92) 
_{May}
(39) 
_{Jun}
(15) 
_{Jul}
(46) 
_{Aug}
(63) 
_{Sep}
(84) 
_{Oct}
(82) 
_{Nov}
(69) 
_{Dec}
(45) 
2016 
_{Jan}
(92) 
_{Feb}
(91) 
_{Mar}
(148) 
_{Apr}
(43) 
_{May}
(58) 
_{Jun}
(117) 
_{Jul}
(92) 
_{Aug}
(140) 
_{Sep}
(49) 
_{Oct}
(33) 
_{Nov}
(85) 
_{Dec}
(40) 
2017 
_{Jan}
(41) 
_{Feb}
(36) 
_{Mar}
(44) 
_{Apr}

_{May}

_{Jun}

_{Jul}

_{Aug}

_{Sep}

_{Oct}

_{Nov}

_{Dec}

S  M  T  W  T  F  S 







1

2

3
(1) 
4

5
(3) 
6

7

8

9

10
(7) 
11
(2) 
12

13

14
(2) 
15
(3) 
16

17

18

19

20

21

22

23
(2) 
24
(2) 
25

26

27

28

29

30

31






From: Roy Stogner <roystgnr@ic...>  20110124 22:38:16

On Mon, 24 Jan 2011, Saurabh Srivastava wrote: >  I'm using Hierarchic basis ...and I checked the second > derivatives for 1st order basis are between 1e5 to 1e11 at several > quadrature points which is though small but nonzero. Indeed on > your suggestion when I changed the basis to Lagrangian i got exact > zeros. Yeah, grep through fe_hierarchic_shape_2D.C for "I have been lazy here". I never needed second derivatives on hierarchic elements for myself, but I didn't want to leave them completely unimplemented, so I just tossed in a central differencing of the first derivatives. That would still come out exactly zero for constant first derivatives, but on your triangles those aren't perfectly constant, because in that case the HIERARCHIC first derivatives are finite differenced too! I'll Cc: this to the mailing lists, in case any altruistic folks want to code up analytic derivative evaluations. Altruistic folks or selfinterested folks  that finite differencing error is probably the limiting factor when we do p refinement on triangles or in 3D.  Roy 
From: yunfei zhu <yzhu2009@go...>  20110124 10:55:16

Dear all, I want to use the ContinuationSystem class. I noted the following comments: enum libMesh::ContinuationSystem::RHS_Mode [protected] There are (so far) two different vectors which may be assembled using the assembly routine: 1.) The Residual = the normal PDE weighted residual 2.) G_Lambda = the derivative wrt the parameter lambda of the PDE weighted residual It is up to the derived class to handle writing separate assembly code for the different cases. Usually something like: switch (rhs_mode) { case Residual: { Fu(i) += ... // normal PDE residual break; } case G_Lambda: { Fu(i) += ... // derivative wrt control parameter break; } My question is that how should I pass the enum rhs_mode to my own assembly function if I write the following code: int main (int argc, char** argv) { //main program } void assemble_jacobian_residual (EquationSystems& es, const std::string& system_name) { assemble the jacobian and residual; //need some code like this switch (rhs_mode) { case Residual: { Fu(i) += ... // normal PDE residual break; } case G_Lambda: { Fu(i) += ... // derivative wrt control parameter break;} } Thank you. Best wishes, yunfei 
From: Roy Stogner <roystgnr@ic...>  20110123 15:28:30

One error jumps right out: for C1 elements you need to provide an analytic gradient function as well.  Roy On Sun, 23 Jan 2011, Ming Q wrote: > I would like to test my simulation on the C1 elements. But when I > try to project the initial values into the system, I got an error. > Could any one tell me how to do the initial values when using a C1 > elements. 
From: Ming Q <micavn@gm...>  20110123 15:10:52

Dear all, I would like to test my simulation on the C1 elements. But when I try to project the initial values into the system, I got an error. Could any one tell me how to do the initial values when using a C1 elements. Thank you in advance, / Ming Q. Here is my code: > //_________________________________________________________________ > Number init_solution (const Point &p, > const Parameters& par, > const std::string& sys_name, > const std::string& unk_name) > { > Real init_value=0.; > Real x = p(0); > Real y = p(1); > if (unk_name.compare("val1")==0) > { > init_value = sin(x)*cos(y); > } > return init_value; > } > //_________________________________________________________________ > void sys1_init (EquationSystems& es, const std::string& sys_name) > { > libmesh_assert ( sys_name == "system1" ); > > TransientLinearImplicitSystem& system=es.get_system<TransientLinearImplicitSystem> ("system1"); > > system.project_solution(init_solution, NULL, es.parameters); > } Then here is the error I got when libMesh doing the projection: > #0 0x00007fff886a7616 in __kill () > #1 0x00007fff88747cca in abort () > #2 0x00000001056e7325 in __gnu_cxx::__verbose_terminate_handler () > #3 0x000000010567bdba in __cxxabiv1::__terminate () > #4 0x000000010567be13 in std::terminate () > #5 0x000000010567be96 in __cxa_throw () > #6 0x000000010065cc9e in libMesh::System::ProjectSolution::operator() (this=0x7fff5fbfe1f0, range=@0x7fff5fbfe1b0) at src/systems/system_projection.C:896 > #7 0x0000000100664d3b in libMesh::Threads::parallel_for<libMesh::StoredRange<libMesh::MeshBase::const_element_iterator, libMesh::Elem const*>, libMesh::System::ProjectSolution> (range=@0x7fff5fbfe1b0, body=@0x7fff5fbfe1f0) at threads.h:270 > #8 0x000000010066348f in libMesh::System::project_vector (this=0x106376e40, fptr=0x100007f6b <init_solution(libMesh::Point const&, libMesh::Parameters const&, std::string const&, std::string const&)>, gptr=0, parameters=@0x7fff5fbfe6c8, new_vector=@0x10639bd00) at src/systems/system_projection.C:282 > #9 0x00000001006637cc in libMesh::System::project_solution (this=0x106376e40, fptr=0x100007f6b <init_solution(libMesh::Point const&, libMesh::Parameters const&, std::string const&, std::string const&)>, gptr=0, parameters=@0x7fff5fbfe6c8) at src/systems/system_projection.C:250 > #10 0x0000000100007f44 in sys1_init (es=@0x7fff5fbfe6c0, sys_name=@0x106376ed8) at source/addinc.cc:56 > #11 0x000000010063c9d1 in libMesh::System::user_initialization (this=0x106376e40) at src/systems/system.C:1505 > #12 0x000000010063cc52 in libMesh::System::init (this=0x106376e40) at src/systems/system.C:204 > #13 0x00000001005d552d in libMesh::EquationSystems::init (this=0x7fff5fbfe6c0) at src/systems/equation_systems.C:119 
From: Roy Stogner <roystgnr@ic...>  20110115 23:58:55

On Sat, 15 Jan 2011, David Fuentes wrote: > for solving a time dependent pde with a spatially varying parameter, > what would be the best data structure to use in libMesh to hold my > spatially varying parameter? > > ie heat equation with thermal conductivity as a function of position, k(x). > > I was thinking to create a separate ExplicitSystem w/ piecewise > constant shape functions across elements. That's pretty much ideal. > However, within a FEMSystem element_time_derivative call, what is the > most efficient way to access the element wise parameters? Not sure. A DiffContext is pretty much restricted to a single system, so you'd have to do something else to access parameters from a separate system. But we don't have an equally easy to use API set up for that access. The obvious solution is for your derived SubFEMSystem to have a pointer to the ExplicitSystem, and then to access that pointer>current_local_solution indexed by the appropriate dof_index... but if you can think of a nicer way to encapsulate that sort of thing, let me know, and we'll see if we can get an API for it into the library.  Roy 
From: David Fuentes <fuentesdt@gm...>  20110115 22:06:50

for solving a time dependent pde with a spatially varying parameter, what would be the best data structure to use in libMesh to hold my spatially varying parameter? ie heat equation with thermal conductivity as a function of position, k(x). I was thinking to create a separate ExplicitSystem w/ piecewise constant shape functions across elements. However, within a FEMSystem element_time_derivative call, what is the most efficient way to access the element wise parameters? thanks, df 
From: Rahul Sampath <rahul.sampath@gm...>  20110115 16:17:55

Hi : I am trying to use LibMesh0.7.0.3 with Petsc3.0.0. I would like to initialize and finalize Petsc and MPI in my own code and not via LibMeshInit. This is mainly because I want to use LibMesh on a subset of communicators of MPI_COMM_WORLD and use Petsc on MPI_COMM_WORLD. I have attached a very simple program that tests this basic functionality. I am encountering a segementation fault or some other MPI_Error when I run this program. I believe this is either a bug in LibMesh or some unsupported feature I am trying to access using this code. Could you kindly look into this issue and give me your suggestions to fix this problem. Thank you. Best regards, Rahul 
From: Roy Stogner <roystgnr@ic...>  20110114 16:17:25

On Fri, 14 Jan 2011, yunfei zhu wrote: > I get warning like this: > WARNING: Second derivatives are not currently correctly calculated on > nonaffine elements! > > I noted that this warning come from function FEBase::compute_map() in > fe_map.C file. > > 00497 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES00498 if > (calculate_d2phi > <http://libmesh.sourceforge.net/doxygen/classlibMesh_1_1FEBase.php#a4170551fa0860822e21475e3c098bc13>)00499 > {00500 libmesh_do_once(00501 libMesh::err > <http://libmesh.sourceforge.net/doxygen/namespacelibMesh.php#ab0817b8542b150bfa3fb23393df99f28>; > << "WARNING: Second derivatives are not currently "00502 > << "correctly calculated on nonaffine elements!"00503 > << std::endl;);00504 }00505 #endif > > > I don't need to calculate the second derivatives for my problem and I did > not enable the second derivatives when I installed libmesh. Yes, you did  they're enabled by default unless you add a configure option to explicitly disable them. If you're not using second derivatives, though, you don't need to worry about the warning. You can try prerequesting only the data you do need from the FE object  if you reinit() the object before requesting data then it'll precompute everything that you might request at a later time, including second derivatives, whether you need it all or not. > But why I get this warning. Because you had an element in your mesh that wasn't detected to be an "affine element" with linear mappings from master space to physical space. That actually shouldn't happen in a mesh that was created from a Tet4>Tet10 conversion unless you moved nodes afterward... if you want to hunt down the nodal coordinates of the offending element and send them to us it would be appreciated.  Roy 
From: yunfei zhu <yzhu2009@go...>  20110114 16:03:41

Dear all, I am working on tracing the equilibrium path of the unstable static nonlinear problem. So, the arclength method are needed. Another scalar parameter \lambda is added to the original equation systems by: system.add_variable("lambda",FIRST,SCALAR); and then solve the augmented system. First, I use element: tet4. It works for small pressure. but when I use mesh.all_second_order() to convert the mesh from tet4 to tet10, I get warning like this: WARNING: Second derivatives are not currently correctly calculated on nonaffine elements! I noted that this warning come from function FEBase::compute_map() in fe_map.C file. 00497 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES00498 if (calculate_d2phi <http://libmesh.sourceforge.net/doxygen/classlibMesh_1_1FEBase.php#a4170551fa0860822e21475e3c098bc13>)00499 {00500 libmesh_do_once(00501 libMesh::err <http://libmesh.sourceforge.net/doxygen/namespacelibMesh.php#ab0817b8542b150bfa3fb23393df99f28>; << "WARNING: Second derivatives are not currently "00502 << "correctly calculated on nonaffine elements!"00503 << std::endl;);00504 }00505 #endif I don't need to calculate the second derivatives for my problem and I did not enable the second derivatives when I installed libmesh. But why I get this warning. By the way, I am trying to using ContinuationSystem to tackle my problem also. I am not totally sure how to use ContinuationSystem until now, but I run some experimental code anyway. This code give the warning above also. I am confused why this happen. could anyone help me? thank you very much. Best wishes, yunfei 
From: Roy Stogner <roystgnr@ic...>  20110111 17:25:02

On Tue, 11 Jan 2011, Vasilis Vavourakis wrote: > load two Meshes: one mesh that is 1D (this is the BEM mesh used for the system) and another dummy one that is 2D (with Trias+Quads+Edges), where > there is onetoone match on all Edge elements. in addition, force the same kind of order of Quadrature for both meshes for surface > integrations...so that there is consistency to quadrature points. this way it can be evaluated the normal vectors...i think! BUT this solution > hinders someone to use AMR directly to the 1D mesh without applying the same thing to the 2D mesh, if you understand what i mean. Actually, if you're using a polygonal/polyhedral domain (or a domain where coarse piecewise quadratic geometry is good enough), the solution here is simple: don't refine the 2D mesh. You'll have to generate "quadrature points" via inverse_map rather than by hoping you matched up quadrature rules, but that's the safe way to go anyway (some of our 2D quadrature rules that you'd need on the boundary of a 3D domain aren't symmetric). A less hackish solution would be to just add (optional) tangent/normal vectors calculation for interior FE evaluations. The FEBase::compute_face_map code is straightforward enough to move and extend; we'd be happy to accept a patch. > Our SparseMatrix class would still parallelize large dense matrices; > you'd be able to get solutions. The only code change is again in > DofMap  the send_list for communication of a BEM variable would need > to be full rather than built up from neighbor topology. > > > right!!! there is no "neighbor topology" in BEM as in FEM ... the DofMap should address all nodes in the mesh!!! > but about DoF index... should the function "dof_number()" do the job? for example: > > Node* node; > ... > const unsigned int index =node>dof_number(system.number(), system.variable_number("u"), 0); > > the above variable "index" is it local or global??? The index is a global index. But getting indices isn't the problem; ensuring sufficient communication between every index is. That's what would require an additional System flag set and DofMap code. > the other problem that arises to this method is that the final system is not of the form: "K * u = F", as in FEM for elasticity; with "K" being > the stiffness and "F" the external force vector. in BEM the final system is like: "H * u = G * t", where matrices "H", "G" contain the boundary > integrals and "u", "t" is the displacement and traction vector to all boundary nodes. so in case there are some boundary nodes to whom the > traction is unknown (known "u" : Dirichlet BC) then there is a need for some sort of system rearrangement!!! which to my mind seems very difficult > for the moment to handle in parallel. This seems to assume that you're solving for u and t in the same function space, and that known coefficients for one correspond to unknowns for the other? In that case you just add_matrix and add_vector to allocate H,G,u_known,t_known terms, assemble them separately, copy the appropriate blocks into your final algebraic system to solve, then copy the results back after the solve.  Roy Requires a little care, but there shouldn't be 
From: Vasilis Vavourakis <vasvav@gm...>  20110111 16:34:52

2011/1/10 Roy Stogner <roystgnr@...> > > On Mon, 10 Jan 2011, Vasilis Vavourakis wrote: > > (1) in BEM mostly boundary integrations occur, which means that in 2D >> problems only 1D discretization is needed. however, since only Edge >> elements >> are "loaded" to the MeshBase class, when i try to evaluate the normal >> vectors to each quadrature point for each boundary (edge) element through >> the "get_normals" function then i get an error. one option would be to >> load >> a 2D mesh (of triangles and quadrilaterals) and perform only boundary >> integration to all nonneighboring side elements, but this way i would >> load >> into the System the nodes located inside the domain. any ideas on that? >> > > The most flexible idea that jumps to my mind: adding a pervariable > flag to identify BEM variables in a System, then querying that in > DofMap so as to only assign degrees of freedom to boundary nodes. > This would be a difficult feature for a nonexperienced libMesh > developer to write, though. > > i've thought of a simple but yet not smart enough solution: load two Meshes: one mesh that is 1D (this is the BEM mesh used for the system) and another dummy one that is 2D (with Trias+Quads+Edges), where there is onetoone match on all Edge elements. in addition, force the same kind of order of Quadrature for both meshes for surface integrations...so that there is consistency to quadrature points. this way it can be evaluated the normal vectors...i think! BUT this solution hinders someone to use AMR directly to the 1D mesh without applying the same thing to the 2D mesh, if you understand what i mean. > > (2) the system matrix "A" of the final linear system "A*x=b" in not banded >> as in FEM, but full populated... this means that parallelization is not >> that >> easy to implement, right? >> > > Our SparseMatrix class would still parallelize large dense matrices; > you'd be able to get solutions. The only code change is again in > DofMap  the send_list for communication of a BEM variable would need > to be full rather than built up from neighbor topology. > right!!! there is no "neighbor topology" in BEM as in FEM ... the DofMap should address all nodes in the mesh!!! but about DoF index... should the function "dof_number()" do the job? for example: Node* node; ... const unsigned int index =node>dof_number(system.number(), system.variable_number("u"), 0); the above variable "index" is it local or global??? > Efficiency is another question, though; I'm not sure what you'd want > to do for preconditioning. > i think that direct solvers will do the job...so i've not come into thinking about preconditioners. the other problem that arises to this method is that the final system is * not* of the form: "K * u = F", as in FEM for elasticity; with "K" being the stiffness and "F" the external force vector. in BEM the final system is like: "H * u = G * t", where matrices "H", "G" contain the boundary integrals and "u", "t" is the displacement and traction vector to all boundary nodes. so in case there are some boundary nodes to whom the traction is unknown (known "u" : Dirichlet BC) then there is a need for some sort of system rearrangement!!! which to my mind seems very difficult for the moment to handle in parallel. Vas >  > Roy > 
From: Roy Stogner <roystgnr@ic...>  20110110 20:57:03

On Tue, 14 Dec 2010, Vetter Roman wrote: > I'm using a 1D mesh with (cubic) Hermite shape functions to > implement a beam. Each element comes with four dof_indices per > variable, namely the nodal value and its derivative for each of the > two nodes. > In beam theory, when using Hermite shape fcts, the derivative dofs > are treated differently than the normal ones in that their stiffness > matrix entries are different. Thus, when building the stiffness > matrix, I need to now which entries of the dof_indices vector > corresponds to which type of dof. > > To be a bit more specific: Currently, I'm iterating over all > active_local_elements. For each of these, I determine the > dof_indices vector, which, for 3 variables, has a size of 12. For > each entry in dof_indices, I would like to know three things: > a) the variable number it belongs to, > b) the node it belongs to, > c) "normal" dof or derivative dof, > so I can build the 12x12 element stiffness matrix correspondingly. > > I'm looking for the "intended" way to do it, i.e. such that my code > would still be correct after doing AMR or even using different shape > functions, both of which I'm planning to do. If you want to be as flexibly correct as possible, the way to do it is to write your formulation in such a way as to avoid needing to know anything about the DoFs "types", just about their shape functions and shape function derivatives. Usually this bites people when they write themselves into a corner with isoparametricLagrangespecific code, but I guess with H2 problems it'd be just as easy to accidentally write cubicHermitespecific code instead. I.e. if you're solving (EI w'')'' = q, then you'd just weight that by a test function and integrate by parts twice, and evaluating the resulting weak formulation should give you a good approximation of w on any C1 function space. > I reckon that the two components per variable (n_comp(s,vn) == 2) > precisely represent the value and the derivative, for each node, > right? (source: > http://www.mailarchive.com/libmeshdevel@.../msg00844.html) Yes. > Is the ordering always guaranteed to be this way? In general, NO. As soon as you bump up to p=4, for example, you'll end up with a new node on which the DoF doesn't correspond to either a value or a derivative. But if you stick to our cubic hermites you should be safe, even with adaptive h refinement.  Roy 
From: Roy Stogner <roystgnr@ic...>  20110110 20:32:41

On Thu, 16 Dec 2010, R P wrote: > I am interested if it is posible to imposse periodicity on the attached mesh (created in Gmsh) in circular (torioidal) direction by this code: > > DofMap & dof_map = system.get_dof_map(); > PeriodicBoundary vert(RealVectorValue(0.0, 0.000001, 0.0)); > > vert.myboundary = 5; > vert.pairedboundary = 6; > > dof_map.add_periodic_boundary(vert); > > If yes, could you please advise me, how to set correct boundary_info IDs (e.g. 5, 6) on the apropriate sides? MeshBase::boundary_info>remove_side() and MeshBase::boundary_info>add_side() are the methods for changing boundary IDs. The best thing you can usually do is loop over all elements, find those which match some criteria (location, face normal, previous boundary ID, etc), and set their new boundary ID accordingly. I've put one of my own little utilities up as an example: http://users.ices.utexas.edu/~roystgnr/meshbcid.C If anyone is interested I could add this utility to the libMesh src/apps/, or accept a patch adding its functionality to src/apps/meshtool.C  Roy 
From: Roy Stogner <roystgnr@ic...>  20110110 20:28:16

On Fri, 17 Dec 2010, robert wrote: > When totally remeshing my domain use a new mesh and EquationSystems > object. If I do this in a loop I would like to copy my EquationSystems > object. > I think about something like: > > time loop{ > > Mesh new_mesh (2); > ...//build some mesh > > EquationSystems new_equation_systems (new_mesh); > ... > ... > old_equation_systems = new_equation_systems; > > > } > > Well, I know that it doesn't work in this way. Is there a > copyconstructor for EquationSystems (I couldn't find it but maybe i > didn't look carefully enough). There isn't  it's such a heavyweight object that we wouldn't want to encourage people to deepcopy it willynilly. > If not, how would you do it? In your case it might be appropriate  if you wanted to write a deep copy constructor (or assignment operator) I'd be happy to look over the patch. > I thought I could write it to a file and read it afterwards again, but > this might not be very efficient. Depending on the I/O type you used, how frequently you remesh, and how expensive your solves in between are, this would either be inefficient or disasterously inefficient.  Roy 
From: Roy Stogner <roystgnr@ic...>  20110110 19:59:18

On Mon, 20 Dec 2010, robert wrote: > I am trying to build a model which solves the same equation as in ex9. > I would like to specify the boundary conditions in a way that the flux > out of the system on the left side equals the flux coming in on the > right side. > I am trying to use something like: > > > for (unsigned int side=0; side<elem>n_sides(); side++) > if (elem>neighbor(side) == NULL) > { > > > const std::vector<std::vector<Real> >& phi_face = fe_face>get_phi(); > > const std::vector<Real>& JxW_face = fe_face>get_JxW(); > > const std::vector<Point >& qface_point = fe_face>get_xyz(); > > fe_face>reinit(elem, side); > > for (unsigned int qp=0; qp<qface.n_points(); qp++) > { > const Real xf = qface_point[qp](0); > const Real yf = qface_point[qp](1); > const Real zf = qface_point[qp](2); > > > > const Real value = neumann_value(xf, yf, zf); > > for (unsigned int i=0; i<phi_face.size(); i++) > Fe(i) += JxW_face[qp]*value*phi_face[i][qp]; > } > > } > but I am not sure how to determine neumann value for the outgoing flux on the other side. Not sure whether this eventually got answered or figured out or not, but basically you want to: Take the xyz value at each quadrature point, Add to it the translation vector between its side and the other side, Call MeshBase::point_locator()() on the result to get the corresponding other_side element. Use FEInterface::inverse_map to get the master element point on the other side, Call fe_other_side>reinit to get flux values there.  Roy 
From: Roy Stogner <roystgnr@ic...>  20110110 19:08:47

On Mon, 3 Jan 2011, Cody Permann wrote: > I have confirmed that there is indeed a problem with sidesets on 3D > meshes when using AMR. Some of the sidesets are completely > incorrect after a single refinement step. The problem can be > reproduced using some of the 3D examples (ex14 or ex15). Any progress on this? If you're seeing it with a modified example can you post a patch to let me reproduce the problem? Thanks,  Roy 
From: Roy Stogner <roystgnr@ic...>  20110110 18:52:46

On Mon, 10 Jan 2011, Vasilis Vavourakis wrote: > (1) in BEM mostly boundary integrations occur, which means that in 2D > problems only 1D discretization is needed. however, since only Edge elements > are "loaded" to the MeshBase class, when i try to evaluate the normal > vectors to each quadrature point for each boundary (edge) element through > the "get_normals" function then i get an error. one option would be to load > a 2D mesh (of triangles and quadrilaterals) and perform only boundary > integration to all nonneighboring side elements, but this way i would load > into the System the nodes located inside the domain. any ideas on that? The most flexible idea that jumps to my mind: adding a pervariable flag to identify BEM variables in a System, then querying that in DofMap so as to only assign degrees of freedom to boundary nodes. This would be a difficult feature for a nonexperienced libMesh developer to write, though. > (2) the system matrix "A" of the final linear system "A*x=b" in not banded > as in FEM, but full populated... this means that parallelization is not that > easy to implement, right? Our SparseMatrix class would still parallelize large dense matrices; you'd be able to get solutions. The only code change is again in DofMap  the send_list for communication of a BEM variable would need to be full rather than built up from neighbor topology. Efficiency is another question, though; I'm not sure what you'd want to do for preconditioning.  Roy 
From: Vasilis Vavourakis <vasvav@gm...>  20110110 17:50:57

hi all and happy new year :) i'd like to ask if anyone has implemented libmesh in Boundary Element Methods. i'm currently trying to built a simple 2D BEM solver for the Laplace PDE "Grad^2(u) = 0". since i'm not that a libmesh expert i face some difficulties in the implementation: (1) in BEM mostly boundary integrations occur, which means that in 2D problems only 1D discretization is needed. however, since only Edge elements are "loaded" to the MeshBase class, when i try to evaluate the normal vectors to each quadrature point for each boundary (edge) element through the "get_normals" function then i get an error. one option would be to load a 2D mesh (of triangles and quadrilaterals) and perform only boundary integration to all nonneighboring side elements, but this way i would load into the System the nodes located inside the domain. any ideas on that? (2) the system matrix "A" of the final linear system "A*x=b" in not banded as in FEM, but full populated... this means that parallelization is not that easy to implement, right? i have more questions but i think that the above are the most "difficult" ones i have so far... in case anyone in the past has worked in that area, their ideas and help would be much appreciated! many thanks, Vas 
From: John Peterson <peterson@ta...>  20110105 20:53:20

On Wed, Jan 5, 2011 at 1:47 PM, Nachiket Gokhale <gokhalen@...> wrote: > Can a quad mesh be generated automatically? I've been told that > triangles and tets can always be generated automatically, while quads > and hexes can't. If you are looking for a proof that a "quadrilateralization" of a given set of points exists (similar to what is known for Delaunay, http://en.wikipedia.org/wiki/Delaunay_triangulation) then I don't think you'll find one. But if you are looking for a practical quadrilateral mesh generation software package, then I think there are a few good ones out there, the OP seems to know of a couple anyway...  John 
From: John Peterson <peterson@ta...>  20110105 19:43:45

On Wed, Jan 5, 2011 at 1:39 PM, Fabrizio De Stefani <deste.fabrizio@...> wrote: > Hi, > I'm trying to build an application which solves an elasticity problem and I > need a mesh generator for the domain discretization. I'd like to know if > libmesh can generate automatically a quad mesh from an arbitrary geometry, > like Salome, Gmhs. Unfortunately, no. Not a quad mesh anyway, just triangles in 2D through an interface to Triangle.  John 
From: Fabrizio De Stefani <deste.fabrizio@gm...>  20110105 19:39:38

Hi, I'm trying to build an application which solves an elasticity problem and I need a mesh generator for the domain discretization. I'd like to know if libmesh can generate automatically a quad mesh from an arbitrary geometry, like Salome, Gmhs. Thanks for your help, Fabrizio 
From: Cody Permann <codypermann@gm...>  20110103 15:00:53

I have confirmed that there is indeed a problem with sidesets on 3D meshes when using AMR. Some of the sidesets are completely incorrect after a single refinement step. The problem can be reproduced using some of the 3D examples (ex14 or ex15). On a regular mesh some of the sidesets will form a "step" pattern cutting into the mesh. Therefore some of the element sides that are actually on the boundary are not in the sideset anymore which is what you are seeing in your problem here. I worked with the problem a little before Christmas after a postdoc reported difficulties to me running some of his 3D problems with AMR and was unable to trace the origin of the problem but I'll take a look at those new lines you mentioned to see if I can run our problems when removing them as you did. What I do know is that the problem has been in there since at least August but I wasn't "easily" able to go back before that point yet. If I learn anything more I'll post to the list again. Cody On Dec 29, 2010, at 3:44 PM, Stefano Berrone wrote: > I recently installed the release 0.7.0.3 (as well as the svn version). I > tried to run a code developed for previous releases and I had some > problems with boundary_id. > The code > > if( elem>neighbor(s) == NULL ) > { > short int side_id = > mesh.boundary_info>boundary_id(elem,s); > ..... > } > > correctly get me the boundary_id set on the boundaries in the initial > mesh up to the first adaptive refinement, but after that it I get > > side_id =1234 > > on a boundary side. > The same code was working with the previous releases. > > I have noticed that some new lines of code are present in the new > version of BoundaryInfo::boundary_id(...). I guess that the problem > comes from that code, in fact if I comment them the code behaves as > expected. > Has anybody had similar problems? > Regards, > Stefano > >  > Learn how Oracle Real Application Clusters (RAC) One Node allows customers > to consolidate database storage, standardize their database environment, and, > should the need arise, upgrade to a full multinode Oracle RAC database > without downtime or disruption > http://p.sf.net/sfu/oraclesfdevnl > _______________________________________________ > Libmeshusers mailing list > Libmeshusers@... > https://lists.sourceforge.net/lists/listinfo/libmeshusers 