project grad f to a VectorFunctionSpace
I don't know if it makes sense but I'd like to do:
m = UnitInterval(11)
Q= FunctionSpace(
W = VectorFunctionS
f = Function(Q)
df = grad(f)
g = project(df,W)
but I get the errpr below. Any help?
Thanks,
Chaffra
/usr/local/
50 Pv = TrialFunction(V)
51 a = ufl.inner(w, Pv)*ufl.dx
---> 52 L = ufl.inner(w, v)*ufl.dx
53
54 # Assemble linear system
/usr/local/
56 if a.shape() == () and b.shape() == ():
57 return a*b
---> 58 return Inner(a, b)
59 #return contraction(a, range(a.rank()), b, range(b.rank()))
60
/usr/local/
160
161 def __new__(cls, a, b):
--> 162 ufl_assert(
163 if isinstance(a, Zero) or isinstance(b, Zero):
164 free_indices, index_dimensions = merge_indices(a, b)
/usr/local/
18 def ufl_assert(
19 "Assert that condition is true and otherwise issue an error with given message."
---> 20 if not condition: error(*message)
21
22
/usr/local/
122 "Write error message and raise an exception."
123 self._log.
--> 124 raise UFLException(
125
126 def begin(self, *message):
UFLException: Shape mismatch.
Question information
- Language:
- English Edit question
- Status:
- Solved
- For:
- DOLFIN Edit question
- Assignee:
- No assignee Edit question
- Solved by:
- Anders Logg
- Solved:
- Last query:
- Last reply:
Revision history for this message
|
#1 |
On Mon, Apr 18, 2011 at 08:10:56PM -0000, Chaffra wrote:
> New question #153352 on DOLFIN:
> https:/
>
> I don't know if it makes sense but I'd like to do:
>
> m = UnitInterval(11)
> Q= FunctionSpace(
> W = VectorFunctionS
>
> f = Function(Q)
> df = grad(f)
>
> g = project(df,W)
>
> but I get the errpr below. Any help?
Skip the dim=3 argument.
--
Anders
> Thanks,
> Chaffra
>
>
> /usr/local/
> 50 Pv = TrialFunction(V)
> 51 a = ufl.inner(w, Pv)*ufl.dx
> ---> 52 L = ufl.inner(w, v)*ufl.dx
> 53
> 54 # Assemble linear system
>
>
> /usr/local/
> 56 if a.shape() == () and b.shape() == ():
> 57 return a*b
> ---> 58 return Inner(a, b)
> 59 #return contraction(a, range(a.rank()), b, range(b.rank()))
>
> 60
>
> /usr/local/
> 160
> 161 def __new__(cls, a, b):
> --> 162 ufl_assert(
> 163 if isinstance(a, Zero) or isinstance(b, Zero):
> 164 free_indices, index_dimensions = merge_indices(a, b)
>
> /usr/local/
> 18 def ufl_assert(
> 19 "Assert that condition is true and otherwise issue an error with given message."
> ---> 20 if not condition: error(*message)
> 21
> 22
>
> /usr/local/
> 122 "Write error message and raise an exception."
> 123 self._log.
> --> 124 raise UFLException(
> 125
> 126 def begin(self, *message):
>
> UFLException: Shape mismatch.
>
>
Revision history for this message
|
#2 |
I tried that and I get this error. The error does not show up for 2D or 3D but for 1D it does. Thanks for you answer. Chaffra
Calling FFC just-in-time (JIT) compiler, this may take some time.
-------
IndexError Traceback (most recent call last)
/home/chaffra/
----> 1
2
3
4
5
/usr/local/
322
323 # Initialize base class
--> 324 MixedFunctionSp
325
326 # FIXME: Add this class:
/usr/local/
298
299 # Initialize base class
--> 300 FunctionSpaceBa
301
302 class VectorFunctionS
/usr/local/
42
43 # JIT-compile element to get ufc_element and ufc_dofmap
---> 44 ufc_element, ufc_dofmap = jit(self.
45
46 # Instantiate DOLFIN FiniteElement and DofMap
/usr/local/
45 # Just call JIT compiler when running in serial
46 if MPI.num_processes() == 1:
---> 47 return local_jit(*args, **kwargs)
48
49 # Compile first on process 0
/usr/local/
122 raise RuntimeError, "Form compiler must implement the jit function."
123
--> 124 return jit_compile(form, parameters=p, common_
125
126
/usr/local/
61 # Check if we get an element or a form
62 if isinstance(object, FiniteElementBase):
---> 63 return jit_element(object, parameters)
64 else:
65 return jit_form(object, parameters, common_cell)
/usr/local/
186
187 # Compile form
--> 188 (compiled_form, module, form_data) = jit_form(form, parameters)
189
190 # Pop cache for element form. Otherwise it might interfere with DOLFIN forms
/usr/local/
131
132 # Generate code
--> 133 compile_
134
135 # Build module using Instant (through UFC)
/usr/local/
138 # Stage 2: intermediate representation
139 cpu_time = time()
--> 140 ir = compute_
141 _print_timing(2, time() - cpu_time)
142
/usr/local/
64 # Compute and flatten representation of integrals
65 info("Computing representation of integrals")
---> 66 irs = [_compute_
67 ir_integrals = [ir for ir in chain(*irs) if not ir is None]
68
/usr/local/
184 form.form_data(),
185 form_id,
--> 186 parameters)
187
188 # Append representation
/usr/local/
57
58 # Compute sum of tensor representations
---> 59 ir["AK"] = _compute_
60
61 elif domain_type == "exterior_facet":
/usr/local/
96 domain_type,
97 facet0, facet1,
---> 98 quadrature_degree)
99
100 # Compute geometry tensor
/usr/local/
26
27 # Compute reference tensor
---> 28 self.A0 = integrate(monomial, domain_type, facet0, facet1, quadrature_order)
29
30 # Extract indices
/usr/local/
48
49 # Compute table Psi for each factor
---> 50 psis = [_compute_psi(v, table, len(points), domain_type) for v in monomial.arguments]
51
52 # Compute product of all Psis
/usr/local/
167 dtuple = _multiindex_
168 # Get values from table
--> 169 Psi[component]
170 else:
171 etable = table[(v.element, v.restriction)]
IndexError: too many indices
Revision history for this message
|
#3 |
On Mon, Apr 18, 2011 at 08:26:56PM -0000, Chaffra wrote:
> Question #153352 on DOLFIN changed:
> https:/
>
> Status: Answered => Open
>
> Chaffra is still having a problem:
> I tried that and I get this error. The error does not show up for 2D or
> 3D but for 1D it does. Thanks for you answer. Chaffra
>
> Calling FFC just-in-time (JIT) compiler, this may take some time.
> -------
> IndexError Traceback (most recent call last)
>
> /home/chaffra/
> ----> 1
> 2
> 3
> 4
> 5
>
> /usr/local/
> 322
> 323 # Initialize base class
>
> --> 324 MixedFunctionSp
> 325
> 326 # FIXME: Add this class:
>
>
> /usr/local/
> 298
> 299 # Initialize base class
>
> --> 300 FunctionSpaceBa
> 301
> 302 class VectorFunctionS
>
> /usr/local/
> 42
> 43 # JIT-compile element to get ufc_element and ufc_dofmap
>
> ---> 44 ufc_element, ufc_dofmap = jit(self.
> 45
> 46 # Instantiate DOLFIN FiniteElement and DofMap
>
>
> /usr/local/
> 45 # Just call JIT compiler when running in serial
>
> 46 if MPI.num_processes() == 1:
> ---> 47 return local_jit(*args, **kwargs)
> 48
> 49 # Compile first on process 0
>
>
> /usr/local/
> 122 raise RuntimeError, "Form compiler must implement the jit function."
> 123
> --> 124 return jit_compile(form, parameters=p, common_
> 125
> 126
>
> /usr/local/
> 61 # Check if we get an element or a form
>
> 62 if isinstance(object, FiniteElementBase):
> ---> 63 return jit_element(object, parameters)
> 64 else:
> 65 return jit_form(object, parameters, common_cell)
>
> /usr/local/
> 186
> 187 # Compile form
>
> --> 188 (compiled_form, module, form_data) = jit_form(form, parameters)
> 189
> 190 # Pop cache for element form. Otherwise it might interfere with DOLFIN forms
>
>
> /usr/local/
> 131
> 132 # Generate code
>
> --> 133 compile_
> 134
> 135 # Build module using Instant (through UFC)
>
>
> /usr/local/
> 138 # Stage 2: intermediate representation
>
> 139 cpu_time = time()
> --> 140 ir = compute_
> 141 _print_timing(2, time() - cpu_time)
> 142
>
> /usr/local/
> 64 # Compute and flatten representation of integrals
>
> 65 info("Computing representation of integrals")
> ---> 66 irs = [_compute_
> 67 ir_integrals = [ir for ir in chain(*irs) if not ir is None]
> 68
>
> /usr/local/
> 184 form.form_data(),
> 185 form_id,
> --> 186 parameters)
> 187
> 188 # Append representation
>
>
> /usr/local/
> 57
> 58 # Compute sum of tensor representations
>
> ---> 59 ir["AK"] = _compute_
> 60
> 61 elif domain_type == "exterior_facet":
>
> /usr/local/
> 96 domain_type,
> 97 facet0, facet1,
> ---> 98 quadrature_degree)
> 99
> 100 # Compute geometry tensor
>
>
> /usr/local/
> 26
> 27 # Compute reference tensor
>
> ---> 28 self.A0 = integrate(monomial, domain_type, facet0, facet1, quadrature_order)
> 29
> 30 # Extract indices
>
>
> /usr/local/
> 48
> 49 # Compute table Psi for each factor
>
> ---> 50 psis = [_compute_psi(v, table, len(points), domain_type) for v in monomial.arguments]
> 51
> 52 # Compute product of all Psis
>
>
> /usr/local/
> 167 dtuple = _multiindex_
> 168 # Get values from table
>
> --> 169 Psi[component]
> 170 else:
> 171 etable = table[(v.element, v.restriction)]
>
> IndexError: too many indices
Could be a bug in FFC. I will have to dig deeper to see what the cause is.
--
Anders
Revision history for this message
|
#4 |
On 18 April 2011 22:10, Chaffra <email address hidden> wrote:
> New question #153352 on DOLFIN:
> https:/
>
> I don't know if it makes sense but I'd like to do:
>
> m = UnitInterval(11)
> Q= FunctionSpace(
> W = VectorFunctionS
>
> f = Function(Q)
> df = grad(f)
In 1D, this will be a scalar, so you don't need the vector space.
Kristian
> g = project(df,W)
>
> but I get the errpr below. Any help?
>
> Thanks,
> Chaffra
>
>
> /usr/local/
> 50 Pv = TrialFunction(V)
> 51 a = ufl.inner(w, Pv)*ufl.dx
> ---> 52 L = ufl.inner(w, v)*ufl.dx
> 53
> 54 # Assemble linear system
>
>
> /usr/local/
> 56 if a.shape() == () and b.shape() == ():
> 57 return a*b
> ---> 58 return Inner(a, b)
> 59 #return contraction(a, range(a.rank()), b, range(b.rank()))
>
> 60
>
> /usr/local/
> 160
> 161 def __new__(cls, a, b):
> --> 162 ufl_assert(
> 163 if isinstance(a, Zero) or isinstance(b, Zero):
> 164 free_indices, index_dimensions = merge_indices(a, b)
>
> /usr/local/
> 18 def ufl_assert(
> 19 "Assert that condition is true and otherwise issue an error with given message."
> ---> 20 if not condition: error(*message)
> 21
> 22
>
> /usr/local/
> 122 "Write error message and raise an exception."
> 123 self._log.
> --> 124 raise UFLException(
> 125
> 126 def begin(self, *message):
>
> UFLException: Shape mismatch.
>
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>
> _______
> Mailing list: https:/
> Post to : <email address hidden>
> Unsubscribe : https:/
> More help : https:/
>
Revision history for this message
|
#5 |
What would be a scalar, df? You should be able to define a vector field on a 1D mesh? It works when you specify dim=2 or 3 but not when dim=1 or dim=None
Revision history for this message
|
#6 |
On 18 April 2011 22:46, Chaffra <email address hidden> wrote:
> Question #153352 on DOLFIN changed:
> https:/
>
> Status: Answered => Open
>
> Chaffra is still having a problem:
> What would be a scalar, df?
Yes, grad(f) == f.dx(0) as far as I know.
>You should be able to define a vector field
> on a 1D mesh?
Yes, that should be possible.
>It works when you specify dim=2 or 3 but not when dim=1
> or dim=None
dim=1 and dim=None will probably be the same since dim will be set
equal to the geometric dimension of the cell, 1D in this case.
I could very well be bug in ffc, I would just expect the MixedElement
from UFL to reduce to a FiniteElement if it contained only one
element.
Kristian
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>
> _______
> Mailing list: https:/
> Post to : <email address hidden>
> Unsubscribe : https:/
> More help : https:/
>
Revision history for this message
|
#7 |
On Mon, Apr 18, 2011 at 08:46:03PM -0000, Chaffra wrote:
> Question #153352 on DOLFIN changed:
> https:/
>
> Status: Answered => Open
>
> Chaffra is still having a problem:
> What would be a scalar, df? You should be able to define a vector field
> on a 1D mesh? It works when you specify dim=2 or 3 but not when dim=1
> or dim=None
Remove the space W and just use Q for the projection.
--
Anders
Revision history for this message
|
#8 |
grad(f) == f.dx(0) only in the dev version of UFL.
Martin
Den 18. apr. 2011 22.58 skrev "Kristian B. Ølgaard" <
<email address hidden>> følgende:
> Question #153352 on DOLFIN changed:
> https:/
>
> Status: Open => Answered
>
> Kristian B. Ølgaard proposed the following answer:
> On 18 April 2011 22:46, Chaffra <email address hidden>
wrote:
>> Question #153352 on DOLFIN changed:
>> https:/
>>
>> Status: Answered => Open
>>
>> Chaffra is still having a problem:
>> What would be a scalar, df?
>
> Yes, grad(f) == f.dx(0) as far as I know.
>
>>You should be able to define a vector field
>> on a 1D mesh?
>
> Yes, that should be possible.
>
>>It works when you specify dim=2 or 3 but not when dim=1
>> or dim=None
>
> dim=1 and dim=None will probably be the same since dim will be set
> equal to the geometric dimension of the cell, 1D in this case.
> I could very well be bug in ffc, I would just expect the MixedElement
> from UFL to reduce to a FiniteElement if it contained only one
> element.
>
> Kristian
>
>>
>> --
>> You received this question notification because you are a member of
>> DOLFIN Team, which is an answer contact for DOLFIN.
>>
>> _______
>> Mailing list: https:/
>> Post to : <email address hidden>
>> Unsubscribe : https:/
>> More help : https:/
>>
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>
> _______
> Mailing list: https:/
> Post to : <email address hidden>
> Unsubscribe : https:/
> More help : https:/
Revision history for this message
|
#9 |
Anders,
Do you mean I should do project(f.dx(0), Q) for the 1D case ( I know that project(grad(f), Q) will not work ) ? I am not very familiar with that notation but if it does the job in 1D that's good for me.
Thanks for all your comments guys.
Revision history for this message
|
#10 |
I think I will treat the 1D case separately. By the way, what is the best way to test for the dimensionality of a mesh, m.topology.dim() or m.geometry().dim()?
Thanks, Chaffra
Revision history for this message
|
#11 |
On Mon, Apr 18, 2011 at 09:33:56PM -0000, Chaffra wrote:
> Question #153352 on DOLFIN changed:
> https:/
>
> Chaffra gave more information on the question:
Projecting to Q works for me, but as Martin points out, this is
something that has been fixed in the development version.
> I think I will treat the 1D case separately. By the way, what is the
> best way to test for the dimensionality of a mesh, m.topology.dim() or
> m.geometry().dim()?
Depends on what you want. The topological dimension is whether you
have intervals, triangles or tetrahedra. The geometric dimension is
the number of coordinates you have for the vertices. The boundary of a
3D tetrahedral mesh is a mesh with topological dimension 2 but
geometric dimension 3 (it's embedded in 3D).
--
Anders
Revision history for this message
|
#13 |
Thanks Anders Logg, that solved my question.