From: Karen L. <kyl...@gm...> - 2010-03-09 21:17:31
|
Dear libmesh users and developers, I have a code for a poisson solver that works on a mesh with about 6000 nodes and about 30k elements. It crashed on the problem that I actually want to solve, which has 1753324 nodes and 11 million elements: ./firstTest -d 3 real.node *** Warning, This code is deprecated, and likely to be removed in future library versions! src/base/libmesh.C, line 349, compiled Feb 25 2010 at 16:11:46 *** found the tetgen files to read Assertion `_id_node.find(foreign_node_id) == _id_node.end()' failed. [0] /Users/karenlee/Documents/Code/libmesh-0.6.4/libmesh/include/mesh/mesh_data.h, line 1042, compiled Feb 25 2010 at 16:14:49 terminate called after throwing an instance of 'libMesh::LogicError' what(): Error in libMesh internal logic Abort trap I tried to go to gdb, said break at main, ran, and got: found the tetgen files to read Assertion `_id_node.find(foreign_node_id) == _id_node.end()' failed. [0] /Users/karenlee/Documents/Code/libmesh-0.6.4/libmesh/include/mesh/mesh_data.h, line 1042, compiled Feb 25 2010 at 16:14:49 terminate called after throwing an instance of 'libMesh::LogicError' what(): Error in libMesh internal logic Program received signal SIGABRT, Aborted. 0x96921e42 in __kill () (gdb) n Single stepping until exit from function __kill, which has no line number information. Program received signal ?, Unknown signal. 0x96921e42 in __kill () (gdb) Quit Any ideas what was wrong and how I can fix it? As far as I can tell, the .node and .ele files are fine... (except they're a bit bit...) Thanks, Karen |
From: Roy S. <roy...@ic...> - 2010-03-09 21:25:42
|
On Tue, 9 Mar 2010, Karen Lee wrote: > Assertion `_id_node.find(foreign_node_id) == _id_node.end()' failed. > [0] > /Users/karenlee/Documents/Code/libmesh-0.6.4/libmesh/include/mesh/mesh_data.h, Looks like void MeshData::add_foreign_node_id() is getting called with the same node id twice. And the code calling it in your case would be TetGenIO::node_in(), which just reads a list of ids and data from your .node file. So either there's a repeated node number somewhere in your file or a bug/incompatibility somewhere in TetGenIO? --- Roy |
From: Karen L. <kyl...@gm...> - 2010-03-09 22:46:13
|
Thanks! Looks like I was just missing some entries of my MeshData since I was combining 2 node files... I'm in gdb mode and it has taken 50 min and it seems that the mesh is still not completely read in yet. Got stalled here: (gdb) n found the tetgen files to read read in nodes and elements n n n n after I pressed n to see if it was actually already done. Is this normal? (1.7M points and 11M elements... Karen On Tue, Mar 9, 2010 at 4:25 PM, Roy Stogner <roy...@ic...>wrote: > > > On Tue, 9 Mar 2010, Karen Lee wrote: > > Assertion `_id_node.find(foreign_node_id) == _id_node.end()' failed. >> [0] >> >> /Users/karenlee/Documents/Code/libmesh-0.6.4/libmesh/include/mesh/mesh_data.h, >> > > Looks like void MeshData::add_foreign_node_id() is getting called with > the same node id twice. And the code calling it in your case would be > TetGenIO::node_in(), which just reads a list of ids and data from your > .node file. So either there's a repeated node number somewhere in > your file or a bug/incompatibility somewhere in TetGenIO? > --- > Roy > |
From: Roy S. <roy...@ic...> - 2010-03-09 22:51:00
|
On Tue, 9 Mar 2010, Karen Lee wrote: > Thanks! Looks like I was just missing some entries of my MeshData since I was combining 2 node files... > > I'm in gdb mode and it has taken 50 min and it seems that the mesh is still not completely read in yet. Got stalled > here: > (gdb) n > found the tetgen files to read > read in nodes and elements > n > n > n > n > > after I pressed n to see if it was actually already done. Is this normal? (1.7M points and 11M elements... Probably. METHOD=dbg turns on a ton of debugging tests, including the GLIBCXX tests, and last time I checked one of those tests turned some basic std::set operations from O(log(N)) to O(N log(N)). This is obviously a problem when N is in the millions... Try METHOD=devel, if you're trying to debug problems that only crop up on large runs. You won't get the vector bounds-checking, which is a shame, but you'll at least still be running all the libmesh_assert() tests, and it'll run much faster. --- Roy |
From: Karen L. <kyl...@gm...> - 2010-03-09 22:52:20
|
I see. Let me try that... I'm looking at example 4, which does not involve reading in a mesh. Is there a way I can parallelize even the read-in? Thanks, Karen On Tue, Mar 9, 2010 at 5:50 PM, Roy Stogner <roy...@ic...>wrote: > > > On Tue, 9 Mar 2010, Karen Lee wrote: > > Thanks! Looks like I was just missing some entries of my MeshData since I >> was combining 2 node files... >> >> I'm in gdb mode and it has taken 50 min and it seems that the mesh is >> still not completely read in yet. Got stalled >> here: >> (gdb) n >> found the tetgen files to read >> read in nodes and elements >> n >> n >> n >> n >> >> after I pressed n to see if it was actually already done. Is this normal? >> (1.7M points and 11M elements... >> > > Probably. METHOD=dbg turns on a ton of debugging tests, including the > GLIBCXX tests, and last time I checked one of those tests turned some > basic std::set operations from O(log(N)) to O(N log(N)). This is > obviously a problem when N is in the millions... > > Try METHOD=devel, if you're trying to debug problems that only crop up > on large runs. You won't get the vector bounds-checking, which is a > shame, but you'll at least still be running all the libmesh_assert() > tests, and it'll run much faster. > --- > Roy > |
From: Roy S. <roy...@ic...> - 2010-03-09 23:07:57
|
On Tue, 9 Mar 2010, Karen Lee wrote: > I see. Let me try that... I'm looking at example 4, which does not > involve reading in a mesh. Is there a way I can parallelize even the > read-in? Yes, but I'm the wrong person to ask; my largest domains are around a million elements, so I haven't yet bothered much with the parallel I/O stuff that Ben and Derek worked on. I think the Nemesis file format is currently our only option for parallel reads/writes. However, even though serial I/O won't scale, in opt mode it shouldn't take much time to begin with. I'd try parallel solves with serial I/O, and if the percentage of your runtime devoted to I/O starts to climb then worry about it then. --- Roy |
From: Karen L. <kyl...@gm...> - 2010-03-09 22:57:41
|
Hm I just did:make METHOD=devel firstTest and I got: Program received signal SIGINT, Interrupt. 0x90ff5980 in __gnu_debug::_Safe_iterator_base::_M_detach () (gdb) quit The program is running. Exit anyway? (y or n) y EIGHT-THREE-NINETY-EIGHT:karen karenlee$ rm firstTest EIGHT-THREE-NINETY-EIGHT:karen karenlee$ make METHOD=devel firstTest echo `../libmesh/contrib/bin/libmesh-config --ldflags` /Users/karenlee/Documents/Code/libmesh-0.6.4/libmesh/lib/i386-apple-darwin9.8.0_devel/libmesh.dylib /Users/karenlee/Documents/Code/libmesh-0.6.4/libmesh/contrib/lib/i386-apple-darwin9.8.0_devel/liblaspack.dylib /Users/karenlee/Documents/Code/libmesh-0.6.4/libmesh/contrib/lib/i386-apple-darwin9.8.0_devel/libmetis.dylib /Users/karenlee/Documents/Code/libmesh-0.6.4/libmesh/contrib/lib/i386-apple-darwin9.8.0_devel/libsfcurves.dylib /Users/karenlee/Documents/Code/libmesh-0.6.4/libmesh/contrib/lib/i386-apple-darwin9.8.0_devel/libgzstream.dylib -lz /Users/karenlee/Documents/Code/libmesh-0.6.4/libmesh/contrib/lib/i386-apple-darwin9.8.0_devel/libgmv.dylib -Wl,-rpath,/opt/local/lib/vtk-5.2 -L/opt/local/lib/vtk-5.2 -lvtkIO -lvtkCommon -lvtkFiltering /Users/karenlee/Documents/Code/libmesh-0.6.4/libmesh/contrib/lib/i386-apple-darwin9.8.0_devel/libtetgen.dylib /Users/karenlee/Documents/Code/libmesh-0.6.4/libmesh/contrib/lib/i386-apple-darwin9.8.0_devel/libtriangle.dylib /Users/karenlee/Documents/Code/libmesh-0.6.4/libmesh/contrib/lib/i386-apple-darwin9.8.0_devel/libexodusii.dylib /Users/karenlee/Documents/Code/libmesh-0.6.4/libmesh/contrib/lib/i386-apple-darwin9.8.0_devel/libnetcdf.dylib /Users/karenlee/Documents/Code/libmesh-0.6.4/libmesh/contrib/lib/i386-apple-darwin9.8.0_devel/libnemesis.dylib /Users/karenlee/Documents/Code/libmesh-0.6.4/libmesh/contrib/lib/i386-apple-darwin9.8.0_devel/libHilbert.dylib -Wl,-rpath,/Users/karenlee/Documents/Code/libmesh-0.6.4/libmesh/lib/i386-apple-darwin9.8.0_devel -Wl,-rpath,/Users/karenlee/Documents/Code/libmesh-0.6.4/libmesh/contrib/lib/i386-apple-darwin9.8.0_devel -L/usr/X11/lib -L/opt/local/lib -L/usr/X11/lib -L/opt/local/lib -Wl,-rpath,/Users/karenlee/Documents/Code/libmesh-0.6.4/libmesh/lib/i386-apple-darwin9.8.0_devel -Wl,-rpath,/Users/karenlee/Documents/Code/libmesh-0.6.4/libmesh/contrib/lib/i386-apple-darwin9.8.0_devel g++ -g -o firstTest `../libmesh/contrib/bin/libmesh-config --include` `../libmesh/contrib/bin/libmesh-config --cxxflags` firstTest.cxx `../libmesh/contrib/bin/libmesh-config --ldflags` i686-apple-darwin9-g++-4.0.1: /Users/karenlee/Documents/Code/libmesh-0.6.4/libmesh/lib/i386-apple-darwin9.8.0_devel/libmesh.dylib: No such file or directory i686-apple-darwin9-g++-4.0.1: /Users/karenlee/Documents/Code/libmesh-0.6.4/libmesh/contrib/lib/i386-apple-darwin9.8.0_devel/liblaspack.dylib: No such file or directory i686-apple-darwin9-g++-4.0.1: /Users/karenlee/Documents/Code/libmesh-0.6.4/libmesh/contrib/lib/i386-apple-darwin9.8.0_devel/libmetis.dylib: No such file or directory i686-apple-darwin9-g++-4.0.1: /Users/karenlee/Documents/Code/libmesh-0.6.4/libmesh/contrib/lib/i386-apple-darwin9.8.0_devel/libsfcurves.dylib: No such file or directory i686-apple-darwin9-g++-4.0.1: /Users/karenlee/Documents/Code/libmesh-0.6.4/libmesh/contrib/lib/i386-apple-darwin9.8.0_devel/libgzstream.dylib: No such file or directory i686-apple-darwin9-g++-4.0.1: /Users/karenlee/Documents/Code/libmesh-0.6.4/libmesh/contrib/lib/i386-apple-darwin9.8.0_devel/libgmv.dylib: No such file or directory i686-apple-darwin9-g++-4.0.1: /Users/karenlee/Documents/Code/libmesh-0.6.4/libmesh/contrib/lib/i386-apple-darwin9.8.0_devel/libtetgen.dylib: No such file or directory i686-apple-darwin9-g++-4.0.1: /Users/karenlee/Documents/Code/libmesh-0.6.4/libmesh/contrib/lib/i386-apple-darwin9.8.0_devel/libtriangle.dylib: No such file or directory i686-apple-darwin9-g++-4.0.1: /Users/karenlee/Documents/Code/libmesh-0.6.4/libmesh/contrib/lib/i386-apple-darwin9.8.0_devel/libexodusii.dylib: No such file or directory i686-apple-darwin9-g++-4.0.1: /Users/karenlee/Documents/Code/libmesh-0.6.4/libmesh/contrib/lib/i386-apple-darwin9.8.0_devel/libnetcdf.dylib: No such file or directory i686-apple-darwin9-g++-4.0.1: /Users/karenlee/Documents/Code/libmesh-0.6.4/libmesh/contrib/lib/i386-apple-darwin9.8.0_devel/libnemesis.dylib: No such file or directory i686-apple-darwin9-g++-4.0.1: /Users/karenlee/Documents/Code/libmesh-0.6.4/libmesh/contrib/lib/i386-apple-darwin9.8.0_devel/libHilbert.dylib: No such file or directory firstTest.cxx: In function ‘void assemble_poisson(EquationSystems&, const std::string&)’: firstTest.cxx:199: warning: unused variable ‘qface_point’ firstTest.cxx:138: warning: unused variable ‘q_point’ make: *** [firstTest] Error 1 The relevant part of my Makefile is as follows: LIBMESH_CONFIG = ../libmesh/contrib/bin/libmesh-config CPPFLAGS = `$(LIBMESH_CONFIG) --include` CXXFLAGS = `$(LIBMESH_CONFIG) --cxxflags` LDFLAGS = LIBS = `$(LIBMESH_CONFIG) --ldflags` PROGS = firstTest all: $(PROGS) firstTest: firstTest.cxx echo $(LIBS) g++ -g -o $@ $(CPPFLAGS) $(CXXFLAGS) firstTest.cxx $(LDFLAGS) $(LIBS) clean: rm -f $(PROGS) Thanks, Karen On Tue, Mar 9, 2010 at 5:52 PM, Karen Lee <kyl...@gm...> wrote: > I see. Let me try that... I'm looking at example 4, which does not involve > reading in a mesh. Is there a way I can parallelize even the read-in? > > Thanks, > Karen > > > On Tue, Mar 9, 2010 at 5:50 PM, Roy Stogner <roy...@ic...>wrote: > >> >> >> On Tue, 9 Mar 2010, Karen Lee wrote: >> >> Thanks! Looks like I was just missing some entries of my MeshData since I >>> was combining 2 node files... >>> >>> I'm in gdb mode and it has taken 50 min and it seems that the mesh is >>> still not completely read in yet. Got stalled >>> here: >>> (gdb) n >>> found the tetgen files to read >>> read in nodes and elements >>> n >>> n >>> n >>> n >>> >>> after I pressed n to see if it was actually already done. Is this normal? >>> (1.7M points and 11M elements... >>> >> >> Probably. METHOD=dbg turns on a ton of debugging tests, including the >> GLIBCXX tests, and last time I checked one of those tests turned some >> basic std::set operations from O(log(N)) to O(N log(N)). This is >> obviously a problem when N is in the millions... >> >> Try METHOD=devel, if you're trying to debug problems that only crop up >> on large runs. You won't get the vector bounds-checking, which is a >> shame, but you'll at least still be running all the libmesh_assert() >> tests, and it'll run much faster. >> --- >> Roy >> > > |
From: Roy S. <roy...@ic...> - 2010-03-09 23:14:45
|
On Tue, 9 Mar 2010, Karen Lee wrote: > Hm I just did:make METHOD=devel firstTest > /Users/karenlee/Documents/Code/libmesh-0.6.4/libmesh/lib/i386-apple-darwin9.8.0_devel/libmesh.dylib: No such file or > directory You've got to do a "METHOD=devel make" in the library itself first. --- Roy |
From: Roy S. <roy...@ic...> - 2010-03-09 23:01:46
|
On Tue, 9 Mar 2010, Karen Lee wrote: > Thanks! Looks like I was just missing some entries of my MeshData since I was combining 2 node files... Oh, but one more suggestion: if you can avoid using MeshData, do so. It doesn't work with more than one MPI rank, it hasn't been tested with threads, and since it hasn't had an active developer in years I'd be wary of it in serial as well. The trouble is that our developers have gotten in the habit of (and in some cases express strong reasons for preferring) using ExplicitSystem solutions to store function data and using subdomain ids to index into tables of discrete data. You seem to be running problems large enough that parallelizing them will be beneficial, and it might be an obstacle if you find yourself having to rewrite the MeshData class before you could run in parallel. --- Roy |
From: Karen L. <kyl...@gm...> - 2010-03-09 23:05:14
|
I don't think I can get away without using MeshData. The reason I'm using libmesh is precisely that I have data on a given (huge) mesh that I am using as input... Without running in parallel then, do you think the problem is too large to converge within a reasonable time frame? Karen On Tue, Mar 9, 2010 at 6:01 PM, Roy Stogner <roy...@ic...>wrote: > > > On Tue, 9 Mar 2010, Karen Lee wrote: > > Thanks! Looks like I was just missing some entries of my MeshData since I >> was combining 2 node files... >> > > Oh, but one more suggestion: if you can avoid using MeshData, do so. > It doesn't work with more than one MPI rank, it hasn't been tested > with threads, and since it hasn't had an active developer in years I'd > be wary of it in serial as well. The trouble is that our developers > have gotten in the habit of (and in some cases express strong reasons > for preferring) using ExplicitSystem solutions to store function data > and using subdomain ids to index into tables of discrete data. You > seem to be running problems large enough that parallelizing them will > be beneficial, and it might be an obstacle if you find yourself having > to rewrite the MeshData class before you could run in parallel. > --- > Roy > |
From: Roy S. <roy...@ic...> - 2010-03-09 23:13:22
|
On Tue, 9 Mar 2010, Karen Lee wrote: > I don't think I can get away without using MeshData. The reason I'm > using libmesh is precisely that I have data on a given (huge) mesh > that I am using as input... Yes, so do many of us, but what kind of data? Data like initial conditions and forcing functions in L2/H1/H2 is easy to store with an ExplicitSystem::solution. Data like material properties maps well to subdomain ids. > Without running in parallel then, do you think the problem is too > large to converge within a reasonable time frame? Impossible to say without knowing what the problem is. If you've got the same code running on a sequence of coarser meshes you're already in a good position to extrapolate runtimes. --- Roy |
From: Karen L. <kyl...@gm...> - 2010-03-10 00:17:10
|
> Yes, so do many of us, but what kind of data? Data like initial > conditions and forcing functions in L2/H1/H2 is easy to store with an > ExplicitSystem::solution. Data like material properties maps well to > subdomain ids. > > It's the gradient of the output of some other simulation, to be used in the RHS of my Poisson's equation. My mesh is the dual mesh of my collaborator who gave me the outputs, and so these numbers are not a smooth function but just whatever comes out. I also just wanted to check how to properly put in the RHS which I read into MeshData. (I have checked that it read in the correct things by outputting it to a file, but not getting the info entry by entry.) Is there a function analogous like get_xyz that I can use rather than doing the following (modified from example 3, and which may not even be correct)? const Real x = q_point[qp](0); const Real y = q_point[qp](1); const Real z = q_point[qp](2); Node node(x,y,z,qp); const Real fxy = mesh_data.operator() (node, 0); Thanks so much, Karen |
From: Roy S. <roy...@ic...> - 2010-03-10 16:38:59
|
On Tue, 9 Mar 2010, Karen Lee wrote: > Yes, so do many of us, but what kind of data? Data like initial > conditions and forcing functions in L2/H1/H2 is easy to store with an > ExplicitSystem::solution. Data like material properties maps well to > subdomain ids. > > It's the gradient of the output of some other simulation, to be used in the RHS of my Poisson's equation. My mesh is > the dual mesh of my collaborator who gave me the outputs, and so these numbers are not a smooth function but just > whatever comes out. If you're not using mesh adaptivity then this is easy, smooth or not - a number associated with every node of your mesh can be stored on a degree 1 LAGRANGE basis, or a number associated with every element can be stored on a degree 0 MONOMIAL basis. If you're doing adaptivity then what you can do depends on what you want associated with new child nodes/elements. > I also just wanted to check how to properly put in the RHS which I read into MeshData. (I have checked that it read > in the correct things by outputting it to a file, but not getting the info entry by entry.) Is there a function > analogous like get_xyz that I can use rather than doing the following (modified from example 3, and which may not > even be correct)? > > const Real x = q_point[qp](0); > const Real y = q_point[qp](1); > const Real z = q_point[qp](2); > > Node node(x,y,z,qp); > > const Real fxy = mesh_data.operator() (node, 0); I think MeshData indexes by nodes, not by points - i.e. if you create a new Node at a quadrature point, or even at the same point as an existing node, MeshData will think you've just handed it a node it doesn't know about. --- Roy |
From: Karen L. <kyl...@gm...> - 2010-03-11 02:14:51
|
> If you're not using mesh adaptivity then this is easy, smooth or not - > a number associated with every node of your mesh can be stored on a > degree 1 LAGRANGE basis, or a number associated with every element can > be stored on a degree 0 MONOMIAL basis. > > If you're doing adaptivity then what you can do depends on what you > want associated with new child nodes/elements. > > > I also just wanted to check how to properly put in the RHS which I read >> into MeshData. (I have checked that it read >> in the correct things by outputting it to a file, but not getting the info >> entry by entry.) Is there a function >> analogous like get_xyz that I can use rather than doing the following >> (modified from example 3, and which may not >> even be correct)? >> >> const Real x = q_point[qp](0); >> const Real y = q_point[qp](1); >> const Real z = q_point[qp](2); >> >> Node node(x,y,z,qp); >> >> >> const Real fxy = mesh_data.operator() (node, 0); >> > > I think MeshData indexes by nodes, not by points - i.e. if you create > a new Node at a quadrature point, or even at the same point as an > existing node, MeshData will think you've just handed it a node it > doesn't know about. > I'm a newbie to finite element, and I guess I wasn't very clear how to store the data on a degree 1 LAGRANGE basis in the code... I just tried to give fxy the correct RHS value and do the same thing as example 3. I know that the quadrature points are not the node values... How do I assign the correct values and not have MeshData think I'm giving it a node it doens't know about? (Also, am I using mesh_data,operator() correctly? And is this the function to use?) Would I be better off storing data to be used in the RHS for elements instead of for the nodes instead? I guess I don't know this too well for now, so might need you to spell it out more... Sorry about that. Thanks! Karen |
From: liang c. <goe...@gm...> - 2010-03-11 05:22:13
|
> > > > > I'm a newbie to finite element, and I guess I wasn't very clear how to > store > the data on a degree 1 LAGRANGE basis in the code... I just tried to give > fxy the correct RHS value and do the same thing as example 3. I know that > the quadrature points are not the node values... How do I assign the > correct > values and not have MeshData think I'm giving it a node it doens't know > about? (Also, am I using mesh_data,operator() correctly? And is this the > function to use?) > > Would I be better off storing data to be used in the RHS for elements > instead of for the nodes instead? > > I guess I don't know this too well for now, so might need you to spell it > out more... Sorry about that. > > Thanks! > Karen > > ------------------------------------------------------------------------------ > > I am newbie also, as far as I know the RHS(right hand side) means the source term in your PDE or weak form of PDE. just my 2cents. Liang |
From: Karen L. <kyl...@gm...> - 2010-03-11 13:42:15
|
Yeah, of course I know that the RHS is my source... At least which equation I wanna solve after all.... I'm just asking how I would write it in the code explicitly, since we no longer have just the vector with indices of the nodes but the points are quadrature points for each element when we construct Fe and Ke... To recap, I guess I wasn't very clear (regarding what Roy mentioned) how to store the data on a degree 1 LAGRANGE basis in the code... I've used the original Poisson example with a constant term for fxy to test the program, but I ought to be feeding that from my data (which I have successfully read into MeshData, just not sure how to access the correct points...) I know that the quadrature points are not the node values... How do I assign the correct values and not have MeshData think I'm giving it a node it doens't know about? Also, should I be using mesh_data.operator()? Would I be better off storing data to be used in the RHS for elements instead of for the nodes instead? Thanks, Karen On Thu, Mar 11, 2010 at 12:22 AM, liang cheng <goe...@gm...> wrote: > >> >> >> I'm a newbie to finite element, and I guess I wasn't very clear how to >> store >> the data on a degree 1 LAGRANGE basis in the code... I just tried to give >> fxy the correct RHS value and do the same thing as example 3. I know that >> the quadrature points are not the node values... How do I assign the >> correct >> values and not have MeshData think I'm giving it a node it doens't know >> about? (Also, am I using mesh_data,operator() correctly? And is this the >> function to use?) >> >> Would I be better off storing data to be used in the RHS for elements >> instead of for the nodes instead? >> >> I guess I don't know this too well for now, so might need you to spell it >> out more... Sorry about that. >> >> Thanks! >> Karen >> >> ------------------------------------------------------------------------------ >> >> I am newbie also, as far as I know the RHS(right hand side) means the > source term in your PDE or weak form of PDE. just my 2cents. > > Liang > |
From: Roy S. <roy...@ic...> - 2010-03-11 15:47:44
|
On Thu, 11 Mar 2010, Karen Lee wrote: > To recap, I guess I wasn't very clear (regarding what Roy mentioned) how to > store the data on a degree 1 LAGRANGE basis in the code... Add a second system, an ExplicitSystem, to your EquationSystems. Add a single variable to it. Load data as a solution into that, and query it when you're integrating your real system. > I've used the original Poisson example with a constant term for fxy > to test the program, but I ought to be feeding that from my data > (which I have successfully read into MeshData, just not sure how to > access the correct points...) I know that the quadrature points are > not the node values... How do I assign the correct values and not > have MeshData think I'm giving it a node it doens't know about? > Also, should I be using mesh_data.operator()? > > Would I be better off storing data to be used in the RHS for elements > instead of for the nodes instead? These are primarily questions about your formulation; figuring out how to implement it is impossible until you've got that pinned down. What is your forcing function, conceptually? A bunch of delta functions at each node? Then you want to integrate your forcing function by looping over all nodes instead of all elements. A smooth or a discontinuous function on your domain? Then you want to integrate over all elements, but interpolate a continuous or a discontinuous finite element solution at each quadrature point. --- Roy |