case 2: return (component == 0 ? p(1) : 0);
// (0, 1-x)
case 3: return (component == 0 ? 0 : 1-p(0));
-
+
// there are
// only four
// shape
template <int dim>
Tensor<1,dim>
FE_Nedelec<dim>::shape_grad_component (const unsigned int i,
- const Point<dim> &/*p*/,
+ const Point<dim> &p,
const unsigned int component) const
{
Assert (i<this->dofs_per_cell, ExcIndexRange(i,0,this->dofs_per_cell));
};
};
+ case 3: // 3d
+ {
+ switch (degree)
+ {
+ // first order Nedelec
+ // elements
+ case 1:
+ {
+ // on the unit cell,
+ // the gradients of
+ // these shape
+ // functions are
+ // linear. we pack them
+ // into an array,
+ // knowing that it may
+ // be expensive to
+ // recompute the whole
+ // array each
+ // time. maybe some
+ // clever compiler can
+ // optimize this out,
+ // seeing that except
+ // for one element all
+ // the other ones are
+ // dead stores...
+ //
+ // the format is: first
+ // index=shape function
+ // number; second
+ // index=vector
+ // component, thrid
+ // index=component
+ // within gradient
+ const double x = p(0),
+ y = p(1),
+ z = p(2);
+ static const double unit_gradients[12][3][3]
+ = { { {0,-(1-z), -(1-y)}, {0,0,0}, { 0, 0, 0} },
+ { {0, 0, 0}, {0,0,0}, { (1-y), -x, 0} },
+ { {0, -z, (1-y)}, {0,0,0}, { 0, 0, 0} },
+ { {0, 0, 0}, {0,0,0}, {-(1-y), -(1-x), 0} },
+
+ { {0, (1-z), -y}, {0,0,0}, { 0, 0, 0} },
+ { {0, 0, 0}, {0,0,0}, { y, x, 0} },
+ { {0, z, y}, {0,0,0}, { 0, 0, 0} },
+ { {0, 0, 0}, {0,0,0}, { -y, (1-x), 0} },
+
+ { {0, 0, 0}, {-(1-z), 0, -(1-x)}, {0, 0, 0} },
+ { {0, 0, 0}, { (1-z), 0, -x}, {0, 0, 0} },
+ { {0, 0, 0}, { z, 0, x}, {0, 0, 0} },
+ { {0, 0, 0}, { -z, 0, (1-x)}, {0, 0, 0} } };
+ // note: simple check
+ // whether this can at
+ // all be: build the
+ // sum over all these
+ // tensors. since the
+ // sum of the shape
+ // functions is a
+ // constant, the
+ // gradient must
+ // necessarily be
+ // zero. this is in
+ // fact the case here,
+ // so test successfull
+ return Tensor<1,dim>(unit_gradients[i][component]);
+ };
+
+ // no other degrees
+ // implemented
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+ };
// presently no other space
// dimension implemented
default:
// for the stokes equation the
// pressure does not couple to
// itself
- std::vector<std::vector<bool> > mask (dim+1, std::vector<bool> (dim+1, true));
- mask[dim][dim] = false;
+// std::vector<std::vector<bool> > mask (dim+1, std::vector<bool> (dim+1, true));
+// mask[dim][dim] = false;
- DoFTools::make_sparsity_pattern (dof_handler, mask, sparsity);
+ DoFTools::make_sparsity_pattern (dof_handler, /*mask,*/ sparsity);
sparsity.compress ();
SparseMatrix<double> A2 (sparsity);
int main ()
{
- std::ofstream logfile("non_primitive_1.output");
+ std::ofstream logfile("non_primitive_2.output");
deallog.attach(logfile);
deallog.depth_console(0);